MATLAB: Help me to solve ordinary differential equation

ode45

Hi,
I am currently working on some codes and I have to find numerical solutions to the following ODE's using Matlab (with graphs):
clc; clear;
mc=0.3;                                                 %g/cm3
wc=0.4;
pc=3.16;                                                %density of cement g/cm3

pw=1;                                                   %density of water g/cm3
T=293;
vc=1/((wc*pc)+1);                                       %the cement volume fraction (dalam bentuk 1 satuan)

ro=(3*vc/(4*pi))^1/3;                                   %jari2 partikel awal (for ex: 2, 4, 8, 16, 32) (dr micron jadi ke mm) tapi masi ratio

L=((4*pi*(wc*pc+1)/3)^(1/3))*ro;                        %nilainya 1   
b1=1231;                                                %K^-1




b2=7579;                                                %K^-1
C=5*10^7;                                               %/mm2h

ER=5364;                                                %K^-1
yg=0.25;                                                %the stoichiometric ratio by mass of water to cement

yw=0.15;                                                %the ratio of water adsorbed in gel pore to reacted cement

RH=0.88;                                                %asumsi

cw=((RH-0.75)/0.25)^3;
v=2;%, ratio of volume of cement gel including gel pore assumed to be 2

rwc3s=0.586;                                            %mass contents of C3S

rwc3a=0.172;                                            %mass contents of C3A

De293=((rwc3s*100)^2.024)*(3.2*10^-14);                 %mm2/h

kr293=(8.05*(10^-10))*((rwc3s*100+rwc3a*100)^0.975);    %mm/h







B293=0.3*10^-10;                                        %mm/h
kr=kr293*exp(-ER*(1/T-1/293));%mm/h
%De=De293*(log(1./alfa)).^1.5;

B=B293*exp(-b1*(1/T-1/293));%mm/h
%kd=(B/alfa^1.5)+C*(Rt-ro).^4; %(mm/h)

tspan = [0 100];   
Rt=0.7;% in the range between 0.5-0.85
syms alfa(t)
input1= diff(alfa,t)==((3*cw)/((yw+yg)*pc*ro^2))*1/((1/(((B/alfa^1.5)+C*(Rt-ro)^4)*ro*alfa^(2/3)))+(((alfa^(-1/3))-((2-alfa)^(-1/3)))/(De293*(log(1/alfa))^1.5))+(1/(kr*ro*alfa^(-2/3))));
cond = alfa(0)==0;
Solve_alfa(t)=ode45(input1,tspan,cond); %or should i use dsolve?
and i also tried to make another code the same but still got NaN value.
clc; clear;
tspan=[0 1000];
y0=[0;0];
 
[t,r]=ode45(@myod,tspan,y0);
plot(t,r);
lgd=legend('r(t)','R(t)');
lgd.FontSize=10;
 
function drdt=myod(t,r) 
wc=0.4;
pc=3.16;                                                %density of cement g/cm3
T=293;
vc=1/((wc*pc)+1);                                       %the cement volume fraction (dalam bentuk 1 satuan)
ro=(3*vc/(4*pi))^1/3;                                   %jari2 partikel awal (for ex: 2, 4, 8, 16, 32) (dr micron jadi ke mm) tapi masi ratio
b1=1231;                                                %K^-1
C=5*10^7;                                               %/mm2h
ER=5364;                                                %K^-1
yg=0.25;                                                %the stoichiometric ratio by mass of water to cement
yw=0.15;                                                %the ratio of water adsorbed in gel pore to reacted cement
RH=0.88;                                                %asumsi
cw=((RH-0.75)/0.25)^3;
v=2;%, ratio of volume of cement gel including gel pore assumed to be 2
rwc3s=0.586;                                            %mass contents of C3S
rwc3a=0.172;                                            %mass contents of C3A
De293=((rwc3s*100)^2.024)*(3.2*10^-14);                 %mm2/h
kr293=(8.05*(10^-10))*((rwc3s*100+rwc3a*100)^0.975);    %mm/h
B293=0.3*10^-10;                                        %mm/h
kr=kr293*exp(-ER*(1/T-1/293));%mm/h
%De=De293*(log(1./alfa)).^1.5;
B=B293*exp(-b1*(1/T-1/293));%mm/h
%kd=(B/alfa^1.5)+C*(Rt-ro).^4; %(mm/h)
Rt=0.75;
    if Rt<=ro
    Sr=0;
    elseif (Rt>=ro)&&(Rt<0.5)
    Sr=4*pi*(Rt)^2;
    elseif (0.5<=Rt)&&(Rt<0.5*(2^0.5))
    Sr=(4*pi*(Rt)^2)-(6*pi*(1-(0.5/(Rt))));
    elseif (Rt>=0.5*(2^0.5)) && (Rt<0.5*(3^0.5))
    syms y x 
    fun = @(y,x) 8*(Rt)./(sqrt((Rt)^2-(x.^2)-(y.^2)));
    ymin=sqrt((Rt)^2-0.5);
    xmin=@(x) sqrt((Rt)^2-0.25-x.^2);
    Sr=integral2(fun,ymin,0.5,xmin,0.5);
    else
    Sr=0;
    end
pw=1;
drdt=zeros(2,1); 
drdt(1)=-((pw*(Sr/(4*pi*(r(2)^2)))*cw)/((yw+yg)*pc*r(1)^2))*1/((1/(((B/(1-(r(1)/ro)^3)^1.5)+C*(r(2)-ro)^4)*r(1)^2))+(((1/r(1))-(1/r(2)))/(De293*(log(1/((1-(r(1)/ro)^3))))^1.5))+(1/(kr*r(1)^2))); 
drdt(2)=(v-1)*(r(1)^2)*(-((pw*(Sr/(4*pi*(r(2)^2)))*cw)/((yw+yg)*pc*r(1)^2))*1/((1/(((B/(1-(r(1)/ro)^3)^1.5)+C*(r(2)-ro)^4)*r(1)^2))+(((1/r(1))-(1/r(2)))/(De293*(log(1/((1-(r(1)/ro)^3))))^1.5))+(1/(kr*r(1)^2))))/Sr;
end
I am new to Matlab and can't seem to figure out how to do these 2 codes. Would it be possible for someone to give me the Matlab codes/instructions I need to find the numerical solutions to these ODE's?
Thank you very much !

Best Answer

The following gets it working, though the output isn't very exciting!. However, you have a number of variables that aren't used.
tspan = [0 100];   
alfa0=0.0001;
[t,alfa]=ode45(@f,tspan,alfa0); %or should i use dsolve?
plot(t,alfa),grid
function dalfadt = f(~,alfa)
%mc=0.3;                         %%%% unused            %g/cm3
wc=0.4;
pc=3.16;                                                %density of cement g/cm3
%pw=1;                           %%%% unused            %density of water g/cm3
T=293;
vc=1/((wc*pc)+1);                                       %the cement volume fraction (dalam bentuk 1 satuan)
ro=(3*vc/(4*pi))^1/3;                                   %jari2 partikel awal (for ex: 2, 4, 8, 16, 32) (dr micron jadi ke mm) tapi masi ratio
%L=((4*pi*(wc*pc+1)/3)^(1/3))*ro;    %%%% unused        %nilainya 1   
b1=1231;                                                %K^-1

%b2=7579;                            %%%% unused        %K^-1
C=5*10^7;                                               %/mm2h
ER=5364;                                                %K^-1
yg=0.25;                                                %the stoichiometric ratio by mass of water to cement
yw=0.15;                                                %the ratio of water adsorbed in gel pore to reacted cement
RH=0.88;                                                %asumsi
cw=((RH-0.75)/0.25)^3;
%v=2;%, %%%% unused  ratio of volume of cement gel including gel pore assumed to be 2
rwc3s=0.586;                                            %mass contents of C3S
rwc3a=0.172;                                            %mass contents of C3A
De293=((rwc3s*100)^2.024)*(3.2*10^-14);                 %mm2/h
kr293=(8.05*(10^-10))*((rwc3s*100+rwc3a*100)^0.975);    %mm/h



B293=0.3*10^-10;                                        %mm/h
kr=kr293*exp(-ER*(1/T-1/293));%mm/h
%De=De293*(log(1./alfa)).^1.5;
B=B293*exp(-b1*(1/T-1/293));%mm/h
%kd=(B/alfa^1.5)+C*(Rt-ro).^4; %(mm/h)
Rt=0.7;% in the range between 0.5-0.85
dalfadt=((3*cw)/((yw+yg)*pc*ro^2))*1/((1/(((B/alfa^1.5)+C*(Rt-ro)^4)*ro*alfa^(2/3)))+(((alfa^(-1/3))-((2-alfa)^(-1/3)))/(De293*(log(1/alfa))^1.5))+(1/(kr*ro*alfa^(-2/3))));
end
Related Question