How to use solve on diff function output

I'm using the diff function to differentiate a symbolic function, but when I try to use the solve function on the resulting expression, I get an incorrect result. This is my code:

syms pmax v az
tau = -pmax/2*((1+v)*(1-az*atan(az^-1))-3/(2*(1+az^2)));
dtau = diff(tau,az);
dfunc = subs(dtau,{az v pmax},{az 0.27 1});
solve(dfunc,az)

which yields this output:

matrix([[5936092008779348386626.6242092545]]).

This is strange to me - I haven't seen this sort of output before and if I plug in a vector and plot the function I get an obvious zero somewhere around 0.5. So what gives?

Here is the code I use for that:

syms pmax v az
tau = -pmax/2*((1+v)*(1-az*atan(az^-1))-3/(2*(1+az^2)));
dtau = diff(tau,az);

x = linspace(0.1,5);

func = subs(tau,{az v pmax},{x 0.27 1});
dfunc = subs(dtau,{az v pmax},{x 0.27 1});

plot(x,func)
hold on
plot(x,dfunc,'r')

I know there might be easier ways to solve this particular problem but this is sort of a systematic problem that I've had over time, and I've done a lot of hacky workarounds, but if anyone knows what I'm doing wrong, could you please help me out?

Thanks!

Answer 1

[Walter Roberson]

dfunc has a 0 at infinity, and 5936092008779348386626.6242092545 is the point where your active Digits setting resulted in a value small enough that the algorithm considered it round-off noise.

One thing to keep in mind is that your 0.27 is only 2 digits of precision, so anything more than 2 digits of precision in the output is numerically unjustified.

You are working with a nonlinear equation for which there is no closed form solution. solve() is defined in that case to attempt a numeric solution, and numeric solutions are defined to only find a single solution if the equation is not a polynomial.

[Steve R]

Thanks for the response.

That means that if I limit the solver over a reasonable range, like [0,10], I can get the zero that I'm looking for, right? Is this possible? I couldn't find anything in the documentation, but it seems like it should be doable.

[Walter Roberson]

Yes, since you know you have a nonlinear equation then there is no point in going for a closed form solution. Instead use vpasolve() and pass in a target range. The solution is at about 0.47092706808875075583019672463