How can i plot dy/dx function?

i have below code. i can get the solution for the derivative equation and able to plot y vs x. But how can i plot dydx vs x?

clear
xspan = (30:0.1:900);
y0 = 0;
[x,y] = ode23(@(x,y) odefun(x,y), xspan, y0);
plot(x,y)

function dydx = odefun(x,y)
R=8.314;

 Q1=0.29*7.54E19*exp(-209201/(R*(x+273)))*(1-y)^3;
 Q2=0.34*5.07E16*exp(-195508/(R*(x+273)))*(1-y)^2;

 
 dydx =(Q1+Q2)/5;
end

Answer 1

[Davide Masiello]

You can use deval in combination with solution structure.

xspan = (30:0.1:900);
y0 = 0;
opts = odeset('AbsTol',1e-8,'RelTol',1e-6);
sol = ode15s(@(x,y) odefun(x,y), xspan, y0, opts);
[~,yp] = deval(sol,sol.x);
plot(sol.x,sol.y)
ylabel('y')
yyaxis right
plot(sol.x,yp)
ylabel('dy/dx')

function dydx = odefun(x,y)
R=8.314;
Q1=0.29*7.54E19*exp(-209201/(R*(x+273)))*(1-y)^3;
Q2=0.34*5.07E16*exp(-195508/(R*(x+273)))*(1-y)^2;

dydx =(Q1+Q2)/5;
end

PS. I changed solver and adjusted tolerances for a better result, you might want to consider doing the same.

Answer 2

[Torsten]

xspan = (30:0.1:900);
y0 = 0;
options = odeset('RelTol',1e-12,'AbsTol',1e-12);
[x,y] = ode15s(@(x,y) odefun(x,y), xspan, y0, options);
dy = zeros(size(y));
for i = 1:numel(x)
  dy(i) = odefun(x(i),y(i));
end
%hold on
%plot(x,y)
plot(x,dy)

%hold off

function dydx = odefun(x,y)
R=8.314;
 Q1=0.29*7.54E19*exp(-209201/(R*(x+273)))*(1-y)^3;
 Q2=0.34*5.07E16*exp(-195508/(R*(x+273)))*(1-y)^2;
 
 dydx =(Q1+Q2)/5;
end

[nur farizan munajat]

Thank you for your answer. Now i expand my function dydx. This funtion is actually based on the experimental data. it suppose to get more than one peak. but i still got only one peak

function dydx = odefun(x,y)
R=8.314;

 Q1=0.29*7.54E19*exp(-209201/(R*(x+273)))*(1-y)^3;
 Q2=0.34*5.07E16*exp(-195508/(R*(x+273)))*(1-y)^2;
 Q3=0.12*1.48E07*exp(-136978/(R*(x+273)))*2*(1-y)*sqrt(-log(1-y));
 Q4=0.09*6.88E07*exp(-142875/(R*(x+273)))*3*(1-y)*(-log(1-y))^(2/3);
 
 dydT =(Q1+Q2+Q3+Q4)/5;
end

[Torsten]

The solution curve for y shows that there should only be one peak ...

Obviously, one reaction dominates the other three massively:

xspan = (30:0.1:900);
y0 = 0;
options = odeset('RelTol',1e-12,'AbsTol',1e-12);
info = 1;
[x1,y1] = ode15s(@(x,y) odefun(x,y,info), xspan, y0, options);
info = 2;
[x2,y2] = ode15s(@(x,y) odefun(x,y,info), xspan, y0, options);
dy1 = zeros(size(y1));
for i = 1:numel(x1)
    dy1(i) = odefun(x1(i),y1(i),1);
end
dy2 = zeros(size(y2));
for i = 1:numel(x2)
    dy2(i) = odefun(x2(i),y2(i),1);
end
figure(1)
hold on
plot(x1,y1)
plot(x2,y2)

Warning: Imaginary parts of complex X and/or Y arguments ignored.

hold off

figure(2)
hold on
plot(x1,dy1)
plot(x2,dy2)

Warning: Imaginary parts of complex X and/or Y arguments ignored.

hold off

function dydT = odefun(x,y, info)
  R=8.314;
  Q1=0.29*7.54E19*exp(-209201/(R*(x+273)))*(1-y)^3;
  Q2=0.34*5.07E16*exp(-195508/(R*(x+273)))*(1-y)^2;
  Q3=0.12*1.48E07*exp(-136978/(R*(x+273)))*2*(1-y)*sqrt(-log(1-y));
  Q4=0.09*6.88E07*exp(-142875/(R*(x+273)))*3*(1-y)*(-log(1-y))^(2/3);
  if info == 1
    dydT = Q1/5;
  else
    dydT = (Q1+Q2+Q3+Q4)/5;
  end
end