Run a system of odes for different parameters and plot

I have a system of 3 odes:

function dydt=odefcnNY_v20(~,y,D,Cs,rho,r0,N,V,k01,F,CL10,V2)
dydt=zeros(3,1);
dydt(1)=(-D*Cs*(1-y(2)))/(rho*r0^2*y(1)); % dr*/dt
dydt(2)=((4*pi*D*N*r0*(1-y(2))*y(1))/V)-(F*k01*y(2)); %dC*/dt
dydt(3)=((k01*y(2)*Cs*V)-(CL10*y(3)))/V2; %dC2/dt
end

I then run this code to simulate the response for a given set of parameters D, Cs, rho, r0, N, V, k01, F, CL10 and V2, to make plots of y(1), y(2), y(3) and various transformations of them and then I compare y(3) to experimental data. Is there a way to run this system of odes in a loop for a different set of k01 values varying from 0.1 to 17000 and to automatically plot the resulting y(3) from the model with the y(3) from the experimental dataset all in one graph, labelled by k01 and clearly labelled by whether it is a model output or an experimental curve?

% Parameters
MW=668.89; % molecular weight 
D=9.916e-5*(MW^-0.4569)*60/600000 %m2/s - [D(cm2/min)=9.916e-5*(MW^-0.4569)*60], divide by 600,000 to convert to m2/s 
rho=1170; %kg/m3 
dv90_1M=3.2e-6; %m 
r0=dv90_1M/2 %m dv90
Cs=1.18e-3 % kg/m3
V=6.9e-8*60;%m3 
W=6.9*1e-6; %kg
N=W/((4/3)*pi*r0^3*rho);
k01=1.7/3600; %s
V2=1.561/1000; % m3
CL10=0.733/3.6e6; % m3/s
F=1;

% initial conditions
tspan=[0 4200*3600]; % 
y0=[1 0 0]; 
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
[t,y]=ode15s(@(t,y) odefcnNY_v20(t,y,D,Cs,rho,r0,N,V,k01,F,CL10,V2), tspan, y0, options);

%Plot y(1) and a transformation of it
plot(t/(3600*24),y(:,1),'-o')%plot time in days, and y(1)
xlabel('time, days')
ylabel('r*, (rp/r0)')
legend('Compound 1')
title ('r*');
plot(t/(3600*24),y(:,1)*r0*1e6); %plot y(1) in microns
xlabel('time, days');
ylabel('r, microns');
legend('Compound 1');
title('r');

%Plot y(2) and a transformation of it
plot(t/(3600*24),y(:,2),'-') %plot time in days, and y(2)
xlabel('time, days')
ylabel('C* (C/Cs)')
legend('Compound 2')
title('C*');
ConcDiss = y(:,2)*Cs; %Concentration dissolved, kg/m3
plot(t/(3600*24), ConcDiss,'DisplayName', 'Theoretical') % time in days, and bulk concentration on y
xlabel('time, days')
ylabel('C1, kg/m3')
title ('C1');

%Plot y(3) and compare to experimental data
plot(t/(3600*24),y(:,3),'-') %plot time in days, and C2
xlabel('time, days')
ylabel('C2')
legend('Compound 1')
title('C2');
hold on
MeasuredData1M=xlsread('PP_Data_1M.xlsx');
t_Measured1M=MeasuredData1M(:,1);
C2_Measured1M=MeasuredData1M(:,2);
plot(t_Measured1M/24,C2_Measured1M,'or','DisplayName','Measured PP1M');
legend('Model', 'Measured PP1M');

[Ameer Hamza]

Can you attach the file 'PP_Data_1M.xlsx'? It will help in understanding how you are trying to visualize this data.

[Nora Rafael]

Sure! Here it is.

Answer 1

[Ameer Hamza]

You can try something like the following code, but there is some issue with your model, which causes ode45 to fail at some point in time.

% Parameters
MW=668.89; % molecular weight 
D=9.916e-5*(MW^-0.4569)*60/600000; %m2/s - [D(cm2/min)=9.916e-5*(MW^-0.4569)*60], divide by 600,000 to convert to m2/s 
rho=1170; %kg/m3 
dv90_1M=3.2e-6; %m 
r0=dv90_1M/2; %m dv90
Cs=1.18e-3; % kg/m3
V=6.9e-8*60;%m3 
W=6.9*1e-6; %kg
N=W/((4/3)*pi*r0^3*rho);
V2=1.561/1000; % m3
CL10=0.733/3.6e6; % m3/s
F=1;
K01 = linspace(0.1, 17000, 5);

MeasuredData1M=xlsread('PP_Data_1M.xlsx');
t_Measured1M=MeasuredData1M(:,1);
C2_Measured1M=MeasuredData1M(:,2);

% initial conditions
tspan=[0 4200*3600]; % 
y0=[1 0 0];

fig = figure();
ax = axes();
hold(ax);

plot(t_Measured1M/24,C2_Measured1M,'or','DisplayName','Measured PP1M');

for i=1:numel(K01)
    k01 = K01(i);
    options = odeset('RelTol',1e-10,'AbsTol',1e-10);
    [t,y]=ode15s(@(t,y) odefcnNY_v20(t,y,D,Cs,rho,r0,N,V,k01,F,CL10,V2), tspan, y0, options);
    plot(t/(3600*24),y(:,3),'-', 'DisplayName', ['model k01=' num2str(k01)]); %plot time in days, and C2
end

%Plot y(3) and compare to experimental data

xlabel('time, days')
ylabel('C2')
title('C2');
legend(ax);

[Nora Rafael]

Thanks!

Indeed something wrong with the code, when r goes to 0 I think.