How to store values from an ode solver?

Hello everyone, I have a huge favor to ask... How in the world can I plot all my values I got from my velocities and reaction rate? The code I have only plots the last value but I need all the values.

function membranereactor
tspan=[0 10] %space time in seconds
f0=[2 0] %molar flow rate fa0 is 2 and fb0 is 0
[t f]=ode45(@membrane,tspan,f0)
global v ra
plot(t,f(:,1),'k-',t,f(:,2),'b-',t,v,'r*',t,ra,'g*','linewidth',1)
legend('Molar Flow Rate of Species A', 'Molar Flow Rate of Species B',...
    'Velocity Profile','Reaction Rate')
xlabel('Space Time (sec)')
title('Membrane Reactor')
grid on
function dFdt=membrane(t,f)
global v ra 
k=0.05 %specific rate constant
kt=0.3 %transport coefficient
ke=0.5 %equilibrium constant
fa0=2 %molar flow rate of A
v0=10
v=v0*((f(1)+f(2))/2) %stoichometric flow change
ca=f(1)/v %concentration of A
cb=f(2)/v %concentration of B
ra=-k*(ca-(cb/ke)) %Reaction rate
dFdt=[ra*v0; (-ra*v0)-((kt*cb)*v0)]

So my end goal is to have 4 curves, the molar flow rates of Fa,Fb,v, and ra as a function of space time. But I cannot get all the values of v and ra to come out. I tried to include the inside the ode but that did not turn out well. I am totally open to new ideas because I am completely stumped.

Thank you for your help, I am super grateful.

Answer 1

[Star Strider]

This works for me:

function dFdt=membrane(t,f)
k=0.05; %specific rate constant
kt=0.3; %transport coefficient
ke=0.5; %equilibrium constant
fa0=2; %molar flow rate of A
v0=10;
v=v0*((f(1)+f(2))/2); %stoichometric flow change
ca=f(1)/v; %concentration of A
cb=f(2)/v; %concentration of B
ra=-k*(ca-(cb/ke)); %Reaction rate
dFdt=[ra*v0; (-ra*v0)-((kt*cb)*v0)];
end

function membranereactor
tspan=[0 10]; %space time in seconds
f0=[2 0]; %molar flow rate fa0 is 2 and fb0 is 0
[t f]=ode45(@membrane,tspan,f0);
k=0.05; %specific rate constant
ke=0.5; %equilibrium constant
v=v0*((f(:,1)+f(:,2))/2); %stoichometric flow change
ca=f(:,1)./v; %concentration of A
cb=f(:,2)./v; %concentration of B
ra=-k*(ca-(cb/ke)); %Reaction rate
figure(1)
plot(t,f(:,1),'k-',t,f(:,2),'b-',t,v,'r*',t,ra,'g*','linewidth',1)
legend('Molar Flow Rate of Species A', 'Molar Flow Rate of Species B',...
    'Velocity Profile','Reaction Rate', 'Location','W')
xlabel('Space Time (sec)')
title('Membrane Reactor')
grid on
end

I eliminated the global calls (using global variables creates code that is difficult to debug so always avoid them), and simply recalculated those values after the ode45 call.

[Marlon Brutus]

Thank you so much. It worked now thanks to you. :D

[Star Strider]

As always, my pleasure!