function f=model2(~,Y,ro)
rt=Y(1);
Rt=Y(2);
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));
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
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;
