fsolve and fimplicit give different results

Hello.

I know there are a lot of questions about nonlinear solver fsolve but I'd like to ask and simultaneously share my experience of the explotatuion of this command.

I want to solve non-linear equation and find a function

where is the digamma function of a complex argument and τ is a real positive number.

Here is my code:

function z=temperature
dd=linspace(0,5,1000);
for j=1:1000
 if j==1
        x0=1;
    else
    x0=temp(j-1);
 end
z1=fsolve(@(y) self(y,dd(j)),x0,optimset('Display','off'));
temp(j)=z1;
end
plot(dd,temp,'Color','black','LineWidth',3);

function z=self(y,x)
tau=1/(10.05);
z=-psin(1/2)+real(psin(1/2+(1+1i)*tau*tanh((1+1i)*x)/y))+log(y);
end
end

Here psin(n,z) for the calculation of the polygamma function for a complex argument can be found here on Matlab central.

The result of this function is the plot

But I can do almost the same with another command fimplicit wit the very simple code

tau=1/(10.05); 
z=@(x,y) -psin(1/2)+real(psin(1/2+(1+1i)*tau*tanh((1+1i)*x)./y))+log(y);
fimplicit(z,[0 5 0 1.1])

and get correct result (see below).

I know about that the solution strongly depends on the initial point for fsolve. But I checked a lot of initial values and I didn't obtain this figure.

Any suggestions?