Plot dB/dt vs. time by ODE45 (with dB/dt on the y axis and time on the x axis)

Hi,

I have plotted the solution for two odes against time.

The two odes are given in the following.

dA/dt=-k1*(exp(-E1/RT))*A

dB/dt=k1*(exp(-E1/RT))*A-k2*(exp(-E2/RT))*B

T=25+5t

E1=60

k1=1e-2

E2=100

k2=1e-2

R=8.314

What I would like to do is also plot dB/dt vs. time (differential equations themselves).

The code I have now is not wokring properly. I was wondering if I make any mistake with the codes?

function diffeqs=odes(t,var)
format long
T=var(1);
A=var(2);
B=var(3);
E1=60;
k1=1e-2;
E2=100;
k2=1e-2;
R=8.314;
 
if t < 20
diffeqs(1,1)=5;%dT/dt
else
diffeqs(1,1)=0;%dT/dt
end
diffeqs(2,1)=-k1*exp(-E1/(R*T))*var(2);%dA/dt
diffeqs(3,1)=k1*exp(-E1/(R*T))*var(2)-k2*exp(-E2/(R*T))*var(3);%dB/dt
end

Plotting

clc;close all; clear all;
timerange=linspace(0,100,50);
ICs= [25,10,0];
[tsol,varsol]=ode45(@(t,var) odes(t,var), timerange, ICs);
pvar = gradient(varsol(:,3));
dt = gradient(tsol);
figure 
yyaxis right
plot(tsol,varsol(:,1),'r','linewidth',2);
xlabel('Time, min');
ylabel('Temperature');
hold on
yyaxis left
plot(tsol,varsol(:,3),'b','linewidth',2);
xlabel('time');
ylabel('B(t)');
grid
legend('B' );
figure
hq = quiver(tsol, varsol(:,3), dt, pvar, 6);
set(hq,'Color','k','linewidth',2)
xlabel('time');
ylabel('dB/dt');
grid
col_header1={'Time','B'}
xlswrite('Book2.xlsx',[tsol(:),varsol(:,3)],'sheet1','A2');
xlswrite('Book2.xlsx',col_header1,'sheet1','A1');

Figure 1 is the plot I got for B(t) vs. time.

Figure 2 is dB/dt vs. time (which is incorrect result)

If I extract the datapoints from Figure 1 and replot it in excel. I got Figure 3 which is the same as Figure 1.

After that, I was using these datapoints to plot dB/dt vs. time then got Figure 4 which is very different from Figure 2.

Therefore, I was wondering if there is anything wrong with my codes? It would be greatly appreciated if anyone can help. Thank you in advance!

Answer 1

[darova]

• What I would like to do is also plot dB/dt vs. time

pvar = gradient(varsol(:,3));
dt = gradient(tsol);
dB = pvar./dt;
plot(tsol,dB)

[S H]

Thanks a lot!