Debugging ode45. Need Help Please.

Here's my code. I'm not sure why it's not running:
function CSTR tspan=[0 10.]; y0=[1.0; 1.0; 100; 80]; % initial concentrations of A, B,Temp Reactor, Temp Cooling Jacket
% A=1.0 M
% B=1.0 M
[t,y]=ode45(@CSTRfunction,tspan,y0);
plot(t,y, 'LineWidth',3);
title('CSTR')
xlabel('Time (mins)'),ylabel('Variables A,B,T, TCJ')
legend('[A]','[B]','T'),grid
end
function dydt=CSTRfunction(time,y)
V= 100;
VR= 100;
E1= -9758.3;
E2= -9758.3;
E3= 8560;
k10= 1.287*10^12;
k20= 1.287*10^12;
k30= 9.043*10^9;
HAB= 4.2;
HBC= -11;
HAD= -41.85;
rho= .9342;
Cp=3.01;
Kw= 4032;
A= .215;
mk= 5;
CpK= 2;
Qk=-1113.5;
k1=k10*exp(E1/(y(3)+273.15));
k2=k20*exp(E2/(y(3)+273.15));
k3=k30*exp(E3/(y(3)+273.15));
dCAdt=V/VR*(y0(1)-y(1))-k1*y(3)*y(1)-k3*y(3)*(y(1))^2;
dCBdt=-V/VR*y(2)+k1*y(3)*y(1)-k2*y(3)*y(2);
dTdt=V/VR*(y0(3)-y(3))-1/rho/Cp*(k1*y(3)*y(1)*HAB+k2*y(3)*y(2)*HBC+k3*y(3)*(y(1))^2*HAD)+Kw*A/rho/Cp/VR*(y(4)-y(3));
dTCJdt=1/mk/CpK*(Qk+Kw*A*(y(3)-y(4)));
% dydt
dydt=[dCAdt,dCBdt,dTRdt,dTRKdt]
end

2 Kommentare

Torsten
Torsten am 29 Jan. 2015
Where do you define y0(1) abd y0(3) ?
Best wishes
Torsten.
Torsten
Torsten am 29 Jan. 2015
In both routines, you will have to define y0 as global.
Best wishes
Torsten.

Melden Sie sich an, um zu kommentieren.

Antworten (0)

Kategorien

Mehr zu General Applications finden Sie in Hilfe-Center und File Exchange

Tags

Gefragt:

am 29 Jan. 2015

Kommentiert:

am 29 Jan. 2015

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by