Plotting output of the ode solver for different initial conditions

Hello,

I am trying to solve set of differential equations. I want change one initial condition and plot all derivaties with respect to the varied initial condition at certain point in timespan. I have tried to run loop for it and added deval into code but this doesn't work. Does someone know how this can be done correctly?

For my example I want to calculate all derivatives at t=0.0000001 and I want to plot all derivates with respect to different T values. At t=0.0000001 for all of them.

Thanks.

Quick summary of my problem:

F vector is a fuction of T and time. I want to take derivative of F with respect to time (dF/dt) and calculate value of F at t=0.0000001 and I want to do same calculation for different T values (500:50:1000). Then, I want to plot F vs T graph. Is it possible ?

for T=500:50:1000 %Set of initial condition
[t,F]=ode15s('che515deneme',[0 0.00001],[1*0.25*2*101.3, 0, 0, 0, 0, 0, 0*0.5*2*101.3, 1*0.75*2*101.3, T, 1])
y=zeros(10,1);
y=deval([t,F],0.0000001);
end

subplot(2,1,1)
plot(T,y(:,2),T,y(:,3),T,y(:,4),T,y(:,5),T,y(:,6),T,y(:,10))
xlabel('T(s)');
ylabel('Coverage')
legend('N*','NH*','NH2*','NH3*','H*','*')

subplot(2,1,2)
plot(t,F(:,1),t,F(:,7),t,F(:,8))
xlabel('t(s)');
ylabel('Coverage')
legend('P_N2','P_NH3','P_H2')

The function that I constructed equations

function dFdt = che515deneme(t,F)
dFdt=zeros(10,1);

%F(1)= P_N2
%F(2)= N*
%F(3)= NH*
%F(4)= NH2*
%F(5)= NH3*
%F(6)= H*
%F(7)= P_NH3
%F(8)= P_H2
%F(9)= T
%F(10)= *


A1f=5.6*10^1;
A1r=2.0*10^10;
A2f=6.0*10^13;
A2r=2.8*10^14;
A3f=4.7*10^13;
A3r=1.8*10^13;
A4f=3.3*10^13;
A4r=9.3*10^12;
A5f=5.9*10^13;
A5r=2.1*10^6;
A6f=5.5*10^5;
A6r=2.3*10^13;

Ea1f=33.0*10^3;
Ea1r=137.0*10^3;
Ea2f=86.5*10^3;
Ea2r=41.2*10^3;
Ea3f=60.4*10^3;
Ea3r=8.6*10^3;
Ea4f=17.2*10^3;
Ea4r=64.6*10^3;
Ea5f=83.7*10^3;
Ea5r=0;
Ea6f=0;
Ea6r=89.4*10^3;

R=8.314;

k1f=A1f*exp(-Ea1f/(R*F(9)));
k1r=A1r*exp(-Ea1r/(R*F(9)));
k2f=A2f*exp(-Ea2f/(R*F(9)));
k2r=A2r*exp(-Ea2r/(R*F(9)));
k3f=A3f*exp(-Ea3f/(R*F(9)));
k3r=A3r*exp(-Ea3r/(R*F(9)));
k4f=A4f*exp(-Ea4f/(R*F(9)));
k4r=A4r*exp(-Ea4r/(R*F(9)));
k5f=A5f*exp(-Ea5f/(R*F(9)));
k5r=A5r*exp(-Ea5r/(R*F(9)));
k6f=A6f*exp(-Ea6f/(R*F(9)));
k6r=A6r*exp(-Ea6r/(R*F(9)));


dFdt(1)=k1r*F(2)^2-k1f*F(1)*F(10)^2;
dFdt(2)=k2r*F(3)*F(10)-k2f*F(2)*F(6)-2*k1r*F(2)^2+2*k1f*F(1)*F(10)^2;
dFdt(3)=k2f*F(2)*F(6)+k3r*F(4)*F(10)-k2r*F(3)*F(10)-k3f*F(3)*F(6);
dFdt(4)=k3f*F(3)*F(6)-k3r*F(4)*F(10)-k4f*F(4)*F(6)+k4r*F(5)*F(10);
dFdt(5)=k4f*F(4)*F(6)-k4r*F(5)*F(10)-k5f*F(5)+k5r*F(7)*F(10);
dFdt(6)=2*k6f*F(8)*F(10)^2-2*k6r*F(6)^2-k2f*F(2)*F(6)+k2r*F(3)*F(10)-k3f*F(3)*F(6)+k3r*F(4)*F(10)-k4f*F(4)*F(6)+k4r*F(5)*F(10);
dFdt(7)=k5f*F(5)-k5r*F(7)*F(10);
dFdt(8)=-k6f*F(8)*F(10)^2+k6r*F(6)^2;
dFdt(10)=-2*k1f*F(1)*F(10)^2+2*k1r*F(2)^2+k2f*F(2)*F(6)-k2r*F(3)*F(10)+k3f*F(3)*F(6)-k3r*F(4)*F(10)+k4f*F(4)*F(6)-k4r*F(5)*F(10)+k5f*F(5)-k5r*F(7)*F(10)-2*k6f*F(8)*F(10)^2+2*k6r*F(6)^2;

Answer 1

[darova]

Here is the solution

[fff fff]

Darova, thank you for help. Now I can collect data from one graph for all initial conditions.

But what I want to get is something different. When I do this, code gives a graph y vs t. I am trying to draw y vs T graph. I want to collect all y (y(1),y(2)...) values at t=0.000005, just a one value not all range between 0 and 0.00001. Then, I want to plot y vs T at t=0.000005. (for example when t=0.000005, y(1)=0.9 at T=500, y(1)=0.8 at T=550, y(1)=0.7 at T=600)

[darova]

I understand. Look

T = 500:50:1000; % Set of initial condition
y = zeros(length(T),10);
for i = 1:length(T)
    [t,F] = ode15s('che515deneme',[0 0.00001],[1*0.25*2*101.3, 0, 0, 0, 0, 0, 0*0.5*2*101.3, 1*0.75*2*101.3, T(i), 1])
    y(i,:) = F(end,:);
end

subplot(2,1,1)
plot(T,y(:,[2:6 10])
xlabel('T(s)');
ylabel('Coverage')
legend('N*','NH*','NH2*','NH3*','H*','*')

subplot(2,1,2)
plot(T,y(:,[1 7 8]))
xlabel('t(s)');
ylabel('Coverage')
legend('P_N2','P_NH3','P_H2')

[fff fff]

Thank you.