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; %kgN=W/((4/3)*pi*r0^3*rho);k01=1.7/3600; %sV2=1.561/1000; % m3CL10=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 onMeasuredData1M=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');
Best Answer