Solving simultaneous ODEs using ode45

I am trying to solve 4 simultaneous ODEs using ode45. One of these coupled ODE is rather complicated, with a heaviside time dependance. I am not very familiar with Matlab so I'm not sure what went wrong when Error messages came up.

The function is:

if true
  % code
end
function output = f33trial(V,F,rho,Cp,Tc0,a,b,T,Fc,Fcs,dFc,T0s,t)
% solving 4 simultaneous ODEs in part 3.3
output=zeros(4,1);
Fc=Fcs+dFc*heaviside(t-5);
output(1)=F/0.7/V*(T0s-T(1))-Fc/0.1/V*(T(1)-Tc0)*(1-exp(-a/rho/Cp*Fc^(b-1)));
output(2)=F/0.1/V*(T(1)-T(2));
output(3)=F/0.1/V*(T(2)-T(3));
output(4)=F/0.1/V*(T(3)-T(4));
end

Here is the script:

if true
  % code
end
% Defining constants
a=5.9e5; % J/(min^0.5dCm^1.5)
b=0.5;
Cp=4187; % J/kgdC
rho=1000; % kg/m^3
F=0.085; % m^3/min
Tc0=25; % dC
V=2.1; % m^3
T0s=150; % dC
Tsps=85; % dC
dFc=0.05; % m^3/min
Fcs=0.51775; % m^3/min obtained from Script32
tspan=[0,10];
Tstart=[T0s,T0s,T0s,T0s];
options=odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
[t,T]=ode45(@(t,T) f33trial(T,t),tspan,Tstart,options);
plot(t,T)

Answer 1

[Torsten]

Try this:

function main
% Defining constants
a=5.9e5; % J/(min^0.5dCm^1.5)
b=0.5;
Cp=4187; % J/kgdC
rho=1000; % kg/m^3
F=0.085; % m^3/min
Tc0=25; % dC
V=2.1; % m^3
T0s=150; % dC
Tsps=85; % dC
dFc=0.05; % m^3/min
Fcs=0.51775; % m^3/min obtained from Script32
tspan=[0,10];
Tstart=[T0s,T0s,T0s,T0s];
options=odeset('RelTol', 1e-6, 'AbsTol', 1e-6);
[T,Y]=ode45(@(t,y) f33trial(t,y,a,b,Cp,rho,F,Tc0,V,T0s,Tsps,dFc,Fcs),tspan,Tstart,options);
plot(T,Y)

function output = f33trial(t,T,a,b,Cp,rho,F,Tc0,V,T0s,Tsps,dFc,Fcs)
output=zeros(4,1);
Fc=Fcs+dFc*heaviside(t-5);
output(1)=F/0.7/V*(T0s-T(1))-Fc/0.1/V*(T(1)-Tc0)*(1-exp(-a/rho/Cp*Fc^(b-1)));
output(2)=F/0.1/V*(T(1)-T(2));
output(3)=F/0.1/V*(T(2)-T(3));
output(4)=F/0.1/V*(T(3)-T(4));
end

Best wishes

Torsten.

[XXC]

Thanks :) It works now but I'm not sure why the original didn't. Also, what is the variable Y (and y) that you put in and why does it make a difference?

[Torsten]

I don't like giving the same name to variables which have different meanings.

I think that the main problem with your code was that your variable lists in your call to f33trial were inconsistent.

Best wishes

Torsten.