It's well known that from time to time nonlinear solver fsolve in Matlab gives controversial results but I'd like to ask and simultaneously share my experience of the exploitation of this command. I want to solve nonlinear equation and find a function $y=y(x)$
$-\Psi(1/2) +{\it \Re} \left( \Psi \left( 1/2+{ \frac { \left( 1+i \right) \tau\,\tanh \left( \left( 1+i \right) x \right) }{y}} \right) \right) +\ln \left( y \right) =0$
where $\Psi(z)$ is the digamma function of a complex argument and $\tau$ is a real positive number. Here is my code in Matlab:
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.
Matlab gives the answer
Meanwhile I can do almost the same with another command fimplicit again in Matlab 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 obtain now correct plot (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. Btw fsolve in Maple solves this nonlinear equation correctly.
How I can improve my code in Matlab with fsolve to reproduce the last figure? Any suggestions?