getting NaN in ODE45

i have a set of nonlinear odes to solve but im getting NaN in output matrix r

i dont know what is wrong with my code, it should work,

ill be more than happy if you share your knowledge on such issue if you've faced before

cheers

here is my function and solver

function dr=Thermal_Bubble2(t,r)

dr=zeros(3,1);

global pa f r0 T0; a=0.150e-6;

l=0.0184;

kb=1.38e-23;

ps=1e5;

c=1500;

d=998;

s=0.0725;

cv=1460;

u=0.000001;

pv=2.33e3;

w=2*pi*f;

b=5.1e-29;

%Equations n=4*pi*(r0^3)*(ps+2*(s/r0)+pv)/(3*kb*T0+(ps+2*(s/r0)+pv)*b) ;

lth=min(r(1)/pi,sqrt(a*r(1)/abs(r(2))));

A=4*pi*(r(1)^3)-3*n*b;

pg=(3*n*kb*r(3))/A;

dr(1)=r(2);

dr(2)=(((1+r(1)/c)/d)*(pg-(2*s/r(1))-(4*u*r(2)/r(1))-ps+pv-pa*sin(w*t))+(r(1)/(d*c))*(((3*n*kb*((4*pi*((r(1))^2)/cv)*((l*(T0-r(3))/lth)-pg*r(2))))/A)-(36*pi*(r(1)^2)*r(2)*n*kb*(r(3))/(A^2))-(2*s*r(2)/(r(1)^2))-(4*u*(((r(2))/(r(2)))^2))-w*pa*cos(w*t))-(1.5*((r(2))^2)*(1-(r(2))/(3*c))))/((r(1))*(1-(r(2))/c));

dr(3)=(4*pi*((r(1))^2)/cv)*((l*(T0-r(3))/lth)-pg*r(2));

solver

clear;

global f pa r0 T0

f=25e3;

pa=50e3;

r0=10e-6;

T0=300;

[t,r]=ode45(@Thermal_Bubble2,(0:0.01/f:20/f),[r0 0 T0]);

plotmatrix(t,r);

Answer 1

[Star Strider]

I believe the problem is here:

dr(2)=(((1+r(1)/c)/d)*(pg-(2*s/r(1))-(4*u*r(2)/r(1))-ps+pv-pa*sin(w*t))+(r(1)/(d*c))*(((3*n*kb*((4*pi*((r(1))^2)/cv)*((l*(T0-r(3))/lth)-pg*r(2))))/A)-(36*pi*(r(1)^2)*r(2)*n*kb*(r(3))/(A^2))-(2*s*r(2)/(r(1)^2))-(4*u*(((r(2))/(r(2)))^2))-w*pa*cos(w*t))-(1.5*((r(2))^2)*(1-(r(2))/(3*c))))/((r(1))*(1-(r(2))/c));
NOTE: r(2)/r(2) in the expression ‘(4*u*(((r(2))/(r(2)))^2))-w*pa*cos(w*t))’

Since ‘r(2)’ is initially zero, this will cause the entire expression to be NaN, which will then propagate throughout your calculations.