How can I solve for this system ODEs?

I am trying to solve for Q,CA,T,Ts for this CSTR. while solving for these system odes and i can't seem to get answers simply using dsolve or ode45. Plz help me.

digits(4)
syms E R T Q(t) CA(t) Ts(t) T(t) t
Qspec=700000; % Specified heat addition
Q0=700000; % heat rate at time 0
CA0=0.25; % concentration of the reactant in the feed
T0=350; % feed temperature
Ts0=350; %
tauTs=10; % time constant for the temperature sensor
tauTH=10; % time constant for actuator dynamics
F=10; %Feed rate
rho=1; % density of the reaction mixture
k0=1.97.*10^24; % rate constant
E/R==20000; %normalized activation energy
Cp=1; % heat capacity of the reaction mixture
delHrxn=160000; %heat of reaction
Vr=100; % volume 
ode1=diff(Q(t),t)==(1./tauTH).*(Qspec-Q(t));
ode2=Vr.*(diff(CA(t),t))==(F.*(CA0-CA(t))./rho)-Vr.*k0.*CA(t).*exp((-20000).*1./T(t));
ode3=Vr.*rho.*Cp.*(diff(T(t),t))==F.*Cp.*(T0-T(t))-Vr.*delHrxn.*CA(t).*k0.*exp((-20000).*1./T(t))+Q(t);
ode4=diff(Ts(t),t)==(1/tauTs).*(T(t)-Ts(t));
eqns=[ode1,ode2,ode3,ode4];
cond1=Q(0)==700000;
cond2=CA(0)==0.25;
cond3=T(0)==350;
cond4=Ts(0)==350;
conds=[cond1 cond2 cond3 cond4];

Using dsolve:

dsolve(eqns,conds)
Warning: Unable to find explicit solution. 
In dsolve (line 201) 
 
ans =
 
[ empty sym ]

Using ode45:

[V] = odeToVectorField(eqns,conds);
M = matlabFunction(V,'vars', {'t','Y'});
ol = ode45(M);


Error using odearguments (line 18)
When the first argument to ode45 is a function handle, the tspan and y0 arguments must be supplied.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

[Walter Roberson]

Your ODE is difficult to solve analytically. In particular, T is difficult.

Your ode45() call needs to pass in a timespan and initial conditions. You should use the two-output form of odeToVectorField in order to see what the order of the variables is so you know the right order to pass in the initial conditions: it will likely be alphabetic, CA Q T Ts .

Numerically, your ode goes near singular in T as you pass t = 8.54

Answer 1

[Star Strider]

Your code has almost done everything you want. You just need to call the ODE solver correctly. (Your system is ‘stiff’, so I chnaged it to ode15s.)

Try this:

syms E R T Q(t) CA(t) Ts(t) T(t) t
Qspec=700000; % Specified heat addition
Q0=700000; % heat rate at time 0
CA0=0.25; % concentration of the reactant in the feed
T0=350; % feed temperature
Ts0=350; %
tauTs=10; % time constant for the temperature sensor
tauTH=10; % time constant for actuator dynamics
F=10; %Feed rate
rho=1; % density of the reaction mixture
k0=1.97.*10^24; % rate constant
E/R==20000; %normalized activation energy
Cp=1; % heat capacity of the reaction mixture
delHrxn=160000; %heat of reaction
Vr=100; % volume 
ode1=diff(Q(t),t)==(1./tauTH).*(Qspec-Q(t));
ode2=Vr.*(diff(CA(t),t))==(F.*(CA0-CA(t))./rho)-Vr.*k0.*CA(t).*exp((-20000).*1./T(t));
ode3=Vr.*rho.*Cp.*(diff(T(t),t))==F.*Cp.*(T0-T(t))-Vr.*delHrxn.*CA(t).*k0.*exp((-20000).*1./T(t))+Q(t);
ode4=diff(Ts(t),t)==(1/tauTs).*(T(t)-Ts(t));
eqns=[ode1,ode2,ode3,ode4];
[V,Subs] = odeToVectorField(eqns);
M = matlabFunction(V,'vars', {'t','Y'});
tspan = [0 50];
ic = [0.25; 7E+4; 350; 350];
[t,y] = ode15s(M, tspan, ic);

lgndc = sym2cell(Subs);
lgndcell = regexp(sprintf('%s\n', lgndc{:}), '\n','split');    % Create Cell Array
figure
plot(t, y)
grid
legend(lgndcell(1:end-1))

Put the ‘conds’ values in the ‘ic’ vector. (The second ‘Subs’ output of odeToVectorField tells the order in which to put them, and also is used in the legend argument.)

EDIT —

Added plot image: