MATLAB: Help me to Solve ODE problem

MATLABmatlab functionode45

Hi every one, I am trying to solve the ode function, but up until now I still couldnot figure it out why my value always NaN.
Here below is my function,
function f=model2(~,Y,ro)
rt=Y(1);
Rt=Y(2);
% Explicit equations
pw=1*0.001;
pc=3.16*0.001;
wc=pw/pc*((1/((0.5^3)*4/3*pi))-1);
T=293;
b1=1231;
C=5*10^7;
ER=5364;
yg=0.25;
yw=0.15;
RH=0.8;
cw=((RH-0.75)/0.25)^3;
v=2;
rwc3s=0.586;
rwc3a=0.172;
De293=((rwc3s*100)^2.024)*(3.2*10^-14);
kr293=(8.05*(10^-10))*((rwc3s*100+rwc3a*100)^0.975);
B293=0.3*10^-10;
alfa=1-(rt/ro)^3;
kr=kr293*exp(-ER*(1/T-1/293));
De=De293*(log(1/alfa))^1.5;
B=B293*exp(-b1*(1/T-1/293));
kd=(B/alfa^1.5)+C*(Rt-ro)^4;
L=((4*pi*(wc*(pc/pw)+1)/3)^(1/3))*ro*0.001;
Rt1=Rt/L;
z=ro/L;
if Rt1<=z
Sr=0;
elseif (Rt1>=z)&&(Rt1<0.5)
Sr=4*pi*(Rt1)^2;
elseif (0.5<=Rt1)&&(Rt1<0.5*(2^0.5))
Sr=(4*pi*(Rt1)^2)-(12*pi*(1-(0.5/(Rt1))));
elseif (Rt1>=0.5*(2^0.5)) && (Rt1<0.5*(3^0.5))
syms y x
fun = @(y,x) 8*(Rt1)/(sqrt((Rt1)^2-(x.^2)-(y.^2)));
ymin=sqrt((Rt1)^2-0.5);
xmin=@(x) sqrt((Rt1)^2-0.25-x.^2);
Sr=integral2(fun,ymin,0.5,xmin,0.5);
else
Sr=0;
end
cst=Sr/(4*pi*(Rt^2));
%Differential equations
drtdt=-((pw*cst*cw)/((yw+yg)*pc*rt^2))*1/((1/(kd*rt^2))+(((1/rt)-(1/Rt))/De)+(1/(kr*rt^2)));
dRtdt=(v-1)*(rt^2)*drtdt/Sr;
f=[drtdt;dRtdt];
end
And here is the command to run the function
clc; clear;
tspan=[0 100];
ro=4*0.001;
y0=[(ro-0.0001) (ro+0.0001)];
[t,y]=ode45(@model2,tspan,y0,[],ro);
figure (2)
plot(t,y(:,1),t,y(:,2));
legend('rt','Rt');
ylabel('mm');
xlabel('t');
axis([tspan 0 100]),grid;
Thank you for you attention! Looking forward to help my question.
Best Regard,
Kevin

Best Answer

Check your equations. Sr is always zero. This causes cst to always be zero, which then causes drdt to also be zero, causing dRdt to be NaN (drdt/Sr = 0/0 = NaN).
Because the output of the first loop is the initial conditions for the next, and you are getting NaN in the first loop, your problem perpetuates. I tried tweaking the initial conditions some, but it didn't help much.
Check why Sr is always 0.