I wonder why NaN appears at Y. How can I fix it?

I was solving this differential equation, but NaN came out. :(

syms Z(R)
g = 9.80665;
%b = ;
r = 0.072;
la = 997;
lb = 1.29;
%k = (la - lb) * g * b^2 / r;
k = 20;
b = sqrt(k * r / (g * (la - lb)));
ode = diff(Z, 2) / (1 + (diff(Z)^2))^1.5 + diff(Z) / (R * (1 + (diff(Z)^2))^2) == 2 + k * Z;
[V] = odeToVectorField(ode);
M = matlabFunction(V, 'vars', {'R', 'Y'});
[R, Y]= ode45(M,[0, 3],[0, 0]);
plot(R, Y);

Answer 1

[Walter Roberson]

Your ode is a function of R, and divides by R, and you started R at 0.

[Walter Roberson]

Look at the second term of the expression. You have defined boundary values at R = 0 with Z'(0)=0. The numerator of the second term is Z' and if we evaluate at R=0 then by assumption in the initial conditions, that is 0 so the numerator would be 0.

Meanwhile the denominator of the second term involves R. At R=0 that would be 0. The other part of the denominator is sqrt(1+Z'^2) and at 0 by way of the initial conditions Z' is 0 there so sqrt(1+0^2) = 1. So the denominator is 0*1

Therefore at R=0 the second term is 0/0 which gives problems. Numeric solvers cannot solve that problem.

Can theory solve the problem? Well, if we understand the equation in terms of limits then if the Z' of the numerator approaches 0 "faster" than the linear R approaches 0 then in the limit the second term is potentially "fast 0 divided by slow 0" which in calculus is 0.

If we temporarily substitute 0 for the second term and look at the first term at 0 with Z'(0)=0 then the first term simplifies to Z''(0) and the right hand side simplifies to 2+k*0 so Z''(0) would have to equal 2. So start exploring equations in which Z''(0)=2, Z' (0)=0, Z(0)=0. Can you make the simplest such equation work with the ode?