Why do i get a empty matrix as a result?

First of all, sorry for the typos/vocabulary, english is not my first.

The thing is, i'm modeling a reactor, with differential equations, that i defined as v(1), v(2) and v(3) for conversion, pressure and temperature. The code run smoothly in general but i need to graph another variable, thats depends on two of these, v(1) and v(2), but the result is a empty matrix.

This is the code

function dv=trioxidoo(t,v)

k1=exp(12.16-5473/v(3)); 
k2=exp(-9.953+8619/v(3)); 
k3=exp(-71.745+52596/v(3)); 
kp=10^(5022/v(3)-4.765); 

pso2=cso20*((thetaso2-v(1))/(1+epsilon*v(1)))*(v(2)/P0)*R*T0;
po2=cso20*((thetao2-0.5*v(1))/(1+epsilon*v(1)))*(v(2)/P0)*R*T0;
pso3=cso20*((thetaso3+v(1))/(1+epsilon*v(1)))*(v(2)/P0)*R*T0;

rso2=-k1*po2*pso2*(1-(pso3/(pso2*(po2^0.5)*kp)))/(22.414*(1+k2*pso2+k3*pso3));

dv1=-rso2/Fso20;
dv2=(-B0/((1-FI)*Ac*Rhoc))*(v(2)/P0)*(v(3)/T0)*(1/(1-epsilon*v(1)));
dv3=((-rso2)*(-Dhrx))/((Fso20*thetaporcp)+(Fso20*v(1)*Dcp));
Ke1=(v(1)/(1-v(1)));
Ke2=(1-(0.5*yso20*v(1)));
Ke3=(yo20-(0.5*yso20*v(1)));
ln=log(Ke1*((Ke2/Ke3)^0.5)*(v(2))^0.5);
Te=(-B/(A+Re*ln));

This are the variables

h=25; %m
D=12; %m

T0=673;% [K]

P0=101.3; % 303.975[kpa]=3atm 1013.25kpa=10atm 1bar=101.3Kpa
R=8.314; %[kpa*m3/kmol*K]

yso20=0.11;
yo20=0.14;
yn20=0.75;
yso30=0;

Mso2=64; %kg/kmol
Mo2=32; %kg/kmol
Mn2=28; %kg/kmol 
Mso3=80; %kg/kmol
PMm=Mso2*yso20+Mo2*yo20+Mn2*yn20+Mso3*yso30; %kg/kmol
rho0=P0*PMm/(R*T0); %kg/m3
v0=350000; %m3/h
m0=v0*rho0; %kg/s
F0=m0/PMm; %kmol/s
Fso20=F0*yso20;

cso20=Fso20/v0; %[kmol/m3]

thetaso2=yso20/yso20;
thetao2=yo20/yso20;
thetan2=yn20/yso20;
thetaso3=yso30/yso20;
delta=(1-0.5-1);
epsilon=delta*yso20;

B0=0.5; %[kpa/m]
FI=0.5;
Ac=(pi/4)*D^2;
Rhoc=3.357;

Dhrx=-99000; %[kJ/kmol]
Cpso2=51.62; %[kJ/kmol*K]
Cpso3=76.87; %[kJ/kmol*K]
Cpo2=33.49; %[kJ/kmol*K]
Cpn2=31.15; %[kJ/kmol*K]
Dcp=Cpso3-0.5*Cpo2-Cpo2; %[kj/kmol*K]
thetaporcp=thetaso2*Cpso2+thetaso3*Cpso3+thetao2*Cpo2+thetan2*Cpn2;%[kj/kmol*K]

wspan=[0:1:50];
[t,v]=ode45('trioxidoo',wspan,[0;P0;T0])

A=0.0938;
B=-98.591;
Re=R;

The variable or result that i need to graph is Te, but as i said earlier, it gives me a empty matrix and i dont know why.

Answer 1

[Sulaymon Eshkabilov]

There are a few errs (where the parameters are defined and how the output variable dv() is defined) in defining the function file trioxidoo.m. Here is the corrected two codes:

% Code to recall and execute trioxidoo.m and plot the computed results: 
clearvars; close all; clc
T0=673;% [K]
P0=101.3; % 303.975[kpa]=3atm 1013.25kpa=10atm 1bar=101.3Kpa
wspan=0:1:50;
[t,v]=ode45('trioxidoo',wspan,[0;P0;T0]);
plot(t,v(:,1), 'r-', t,v(:,2), 'g-', t,v(:,3), 'b-','linewidth', 2)
grid on; xlabel('t'); ylabel('v_1, v_2, v_3')
legend('v_1', 'v_2', 'v_3', 'location', 'best')
%% You can also consider to plot the solutions in this way.
% Since the magnitude of v_1 w.r.t v_2 and v_3 is so small.
figure
yyaxis left
plot(t,v(:,1), 'b-', 'linewidth', 2)
ylabel('v_1')
grid on
yyaxis right
plot(t,v(:,2), 'r--', t,v(:,3), 'g-.','linewidth', 2)
ylabel('v_2, v_3'), ylim([50, 1000])
xlabel('t')
legend('v_1', 'v_2', 'v_3', 'location', 'best')

Here is the corrected trioxidoo.m function file:

function dv=trioxidoo(t,v)
h=25; %m
D=12; %m

T0=673;% [K]

P0=101.3; % 303.975[kpa]=3atm 1013.25kpa=10atm 1bar=101.3Kpa
R=8.314; %[kpa*m3/kmol*K]

yso20=0.11;
yo20=0.14;
yn20=0.75;
yso30=0;

Mso2=64; %kg/kmol
Mo2=32; %kg/kmol
Mn2=28; %kg/kmol 
Mso3=80; %kg/kmol
PMm=Mso2*yso20+Mo2*yo20+Mn2*yn20+Mso3*yso30; %kg/kmol
rho0=P0*PMm/(R*T0); %kg/m3
v0=350000; %m3/h
m0=v0*rho0; %kg/s
F0=m0/PMm; %kmol/s
Fso20=F0*yso20;

cso20=Fso20/v0; %[kmol/m3]

thetaso2=yso20/yso20;
thetao2=yo20/yso20;
thetan2=yn20/yso20;
thetaso3=yso30/yso20;
delta=(1-0.5-1);
epsilon=delta*yso20;

B0=0.5; %[kpa/m]
FI=0.5;
Ac=(pi/4)*D^2;
Rhoc=3.357;

Dhrx=-99000; %[kJ/kmol]
Cpso2=51.62; %[kJ/kmol*K]
Cpso3=76.87; %[kJ/kmol*K]
Cpo2=33.49; %[kJ/kmol*K]
Cpn2=31.15; %[kJ/kmol*K]
Dcp=Cpso3-0.5*Cpo2-Cpo2; %[kj/kmol*K]
thetaporcp=thetaso2*Cpso2+thetaso3*Cpso3+thetao2*Cpo2+thetan2*Cpn2;%[kj/kmol*K]
A=0.0938;
B=-98.591;
Re=R;
k1=exp(12.16-5473/v(3)); 
k2=exp(-9.953+8619/v(3)); 
k3=exp(-71.745+52596/v(3)); 
kp=10^(5022/v(3)-4.765); 

pso2=cso20*((thetaso2-v(1))/(1+epsilon*v(1)))*(v(2)/P0)*R*T0;
po2=cso20*((thetao2-0.5*v(1))/(1+epsilon*v(1)))*(v(2)/P0)*R*T0;
pso3=cso20*((thetaso3+v(1))/(1+epsilon*v(1)))*(v(2)/P0)*R*T0;

rso2=-k1*po2*pso2*(1-(pso3/(pso2*(po2^0.5)*kp)))/(22.414*(1+k2*pso2+k3*pso3));

dv(1,:)=-rso2/Fso20;
dv(2,:)=(-B0/((1-FI)*Ac*Rhoc))*(v(2)/P0)*(v(3)/T0)*(1/(1-epsilon*v(1)));
dv(3,:)=((-rso2)*(-Dhrx))/((Fso20*thetaporcp)+(Fso20*v(1)*Dcp));
Ke1=(v(1)/(1-v(1)));
Ke2=(1-(0.5*yso20*v(1)));
Ke3=(yo20-(0.5*yso20*v(1)));
ln=log(Ke1*((Ke2/Ke3)^0.5)*(v(2))^0.5);
Te=(-B/(A+Re*ln));
end