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 Comments

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

Sign in to comment.

Answers (0)

Categories

Find more on General Applications in Help Center and File Exchange

Tags

Asked:

on 29 Jan 2015

Commented:

on 29 Jan 2015

Community Treasure Hunt

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

Start Hunting!