plotting profiles from differential equations

Can anyone help me correct my code to produce a graphical solution to my differrential equations? Below is my code which is saved as Reactor3.m. I have attached my differretial equations in a pdf file along with the intial conditions.

function [dFdV,dTdV] = Reactor3(V,F,T)
%Reactor3 function is used to obtain mass and energy balances across the
%reactor 

%Constants
P=2200; %inlet Pressure (kPa)
dcat=850000; %catalyst density g/m3
U=0.17603*3600; %overall heat transfer co (kj/m2 h K)
a=80; %(m^-1)
Tc=433.15; %coolant temperature isothermal (K)
CpA=64.7081297; %kj/kmol k
CpB=33.43828781; %kj/kmol k
CpC=73.79376; %kj/kmol k
CpD=45.25448885; %kj/kmol k
CpE=35.69429546; %kj/kmol k
Cp_arg= 21.08524537; %kj/kmol k
Cp_meth=48.47777972; %kj/kmol k
Cp_eth=81.59024981; %kj/kmol k
Cp_N2=29.91277101; %kj/kmol k
F_arg=2545.557663; %kmol/h
F_meth=25855.7708; %kmol/h
F_eth=342.8308447; %kmol/h
F_N2=3011.063741; %kmol/h
F_DCE=0.009610; %kmol/h
H1=-117.152*1000; %kj/kmol
H2=-1322.144*1000; %kj/kmol

%Kinetic Constants
k1=6.867*exp(-9260/(1.987*T(1)));
k2=107.3*exp(-10413/(1.987*T(1)));
k3=10.62*exp(-8792/(1.987*T(1)));
k4=39.59*exp(-9473/(1.987*T(1)));
k5=0.000001346*exp(9610/(1.987*T(1)));
k6=0.00000003109*exp(15044/(1.987*T(1)));
k7=0.19;
k8=0.07;

Ft = F(1)+F(2)+F(3)+F(4)+F(5)+ F_arg+ F_meth+ F_N2+ F_DCE;  %total molar flowrate

P_O2 = P*( F(2)/Ft); %partial pressure of O2 (kPa)
P_E = P*( F(1)/Ft); %partial pressure of ethylene (kPa)
P_DCE= P*(F_DCE/Ft); %partial pressure of DCE (kPa)

r1=(dcat*((k1*P_O2*P_E)-(k2*P_O2*P_E*P_DCE^k7)))/1000*(1+k5*P_O2+k6*P_E); %kmol/ h m3
r2=(dcat*((k3*P_O2*P_E)-(k4*P_O2*P_E*P_DCE^k8)))/1000*(1+k5*P_O2+k6*P_E); %kmol/ h m3

%flowrate * heat capacity of each component 
C_A= F(1)*CpA; %kj/h K
C_B= F(2)*CpB; %kj/h K
C_C= F(3)*CpC; %kj/h K
C_D= F(4)*CpD; %kj/h K
C_E= F(5)*CpE; %kj/h K
C_meth=F_meth*Cp_meth; %kj/h K
C_eth=F_eth*Cp_eth; %kj/h K
C_N2=F_N2*Cp_N2; %kj/h K
C_arg=F_arg*Cp_arg; %kj/h K

%sum of heat capacity flowrates
FC= C_A+C_B+C_C+C_D+C_E+C_meth+C_eth+C_N2+C_arg; %kj/h K

%Mass balance
dFdV(1) = r1+r2;
dFdV(2) = 0.5*r1+3*r2;
dFdV(3) = -r1;
dFdV(4) = -2*r2;
dFdV(5) = -2*r2;

%Energy balance
dTdV(1) = (U*a*(T-Tc)+dcat*(H1*r1+H2*r2))./FC;
end

[V,T] = ode45('Reactor3', [0 100], [523.15]);

plot(V,T)
xlabel('V(m3)');
ylabel('T (K))');

[Alex Mcaulley]

What values of V, F and T are you using? What is the error are you getting?

[Liam Gavin]

My values of V range from 0-100 m3, and im trying to deduce how the flowrate (F) and teperature (T) change across my volume.

Errors i get when i run the reactor3 file are:

Not enough input arguments.

Error in Reactor3 (line 29)

k1=6.867*exp(-9260/(1.987*T(1)));

Whereas the errors i get when i try to plot the graph include:

Index exceeds array bounds.

Error in Reactor3 (line 29)

k1=6.867*exp(-9260/(1.987*T(1)));

Error in odearguments (line 90)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)

odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in energy_balance_solver (line 1)

[V,T] = ode45('Reactor3', [0 100], [523.15]);

Answer 1

[Torsten]

function dFTdV = Reactor3(V,FT)

%Reactor3 function is used to obtain mass and energy balances across the
%reactor 
F = FT(1:5);
T(1) = FT(6);

%Constants
P=2200; %inlet Pressure (kPa)
dcat=850000; %catalyst density g/m3
U=0.17603*3600; %overall heat transfer co (kj/m2 h K)
a=80; %(m^-1)
Tc=433.15; %coolant temperature isothermal (K)
CpA=64.7081297; %kj/kmol k
CpB=33.43828781; %kj/kmol k
CpC=73.79376; %kj/kmol k
CpD=45.25448885; %kj/kmol k
CpE=35.69429546; %kj/kmol k
Cp_arg= 21.08524537; %kj/kmol k
Cp_meth=48.47777972; %kj/kmol k
Cp_eth=81.59024981; %kj/kmol k
Cp_N2=29.91277101; %kj/kmol k
F_arg=2545.557663; %kmol/h
F_meth=25855.7708; %kmol/h
F_eth=342.8308447; %kmol/h
F_N2=3011.063741; %kmol/h
F_DCE=0.009610; %kmol/h
H1=-117.152*1000; %kj/kmol
H2=-1322.144*1000; %kj/kmol

%Kinetic Constants
k1=6.867*exp(-9260/(1.987*T(1)));
k2=107.3*exp(-10413/(1.987*T(1)));
k3=10.62*exp(-8792/(1.987*T(1)));
k4=39.59*exp(-9473/(1.987*T(1)));
k5=0.000001346*exp(9610/(1.987*T(1)));
k6=0.00000003109*exp(15044/(1.987*T(1)));
k7=0.19;
k8=0.07;
Ft = F(1)+F(2)+F(3)+F(4)+F(5)+ F_arg+ F_meth+ F_N2+ F_DCE;  %total molar flowrate
P_O2 = P*( F(2)/Ft); %partial pressure of O2 (kPa)
P_E = P*( F(1)/Ft); %partial pressure of ethylene (kPa)
P_DCE= P*(F_DCE/Ft); %partial pressure of DCE (kPa)
r1=(dcat*((k1*P_O2*P_E)-(k2*P_O2*P_E*P_DCE^k7)))/1000*(1+k5*P_O2+k6*P_E); %kmol/ h m3
r2=(dcat*((k3*P_O2*P_E)-(k4*P_O2*P_E*P_DCE^k8)))/1000*(1+k5*P_O2+k6*P_E); %kmol/ h m3

%flowrate * heat capacity of each component 
C_A= F(1)*CpA; %kj/h K
C_B= F(2)*CpB; %kj/h K
C_C= F(3)*CpC; %kj/h K
C_D= F(4)*CpD; %kj/h K
C_E= F(5)*CpE; %kj/h K
C_meth=F_meth*Cp_meth; %kj/h K
C_eth=F_eth*Cp_eth; %kj/h K
C_N2=F_N2*Cp_N2; %kj/h K
C_arg=F_arg*Cp_arg; %kj/h K

%sum of heat capacity flowrates
FC= C_A+C_B+C_C+C_D+C_E+C_meth+C_eth+C_N2+C_arg; %kj/h K

%Mass balance
dFdV(1) = r1+r2;
dFdV(2) = 0.5*r1+3*r2;
dFdV(3) = -r1;
dFdV(4) = -2*r2;
dFdV(5) = -2*r2;

%Energy balance
dTdV(1) = (U*a*(T-Tc)+dcat*(H1*r1+H2*r2))./FC;

dFTdV = [dFdV.';dTdV];
end

[V,FT] = ode15s(@Reactor3, [0 100], [13096.019;3156.3;0;32.38;0;523.15]);

plot(V,FT(:,6))
xlabel('V(m3)');
ylabel('T (K))');

[Liam Gavin]

Hi, I really appreciate the help!

However the variables F and T are differrent. I know that all differential equations are dependant on T and most of them are dependant on F.

Is there anyway that i could plot how just T changes with V? Not FT?

[Torsten]

T depends on F ; thus you have to solve all 6 equations simultaneoulsy.

The plot shows how T changes with V. I don't understand what you mean by your last question.