MATLAB: ODE Solver wont input P1 function

ode45

function first_oder_ode 
 
% SOLVE  dx/dt = -3 exp(-t). 
% initial conditions: x(0) = 0 
CV1=110*0.1336806;
CV2=110*0.1336806;
Area=15; % cross sectional area of tank, m^2
P0=1.0*10^5;
P3=1.1*10^5;
rho=1000;
g=9.8;
k0=10.09782398;
kM=423463.7642;
ka=8.869925491;
k1=0.633075579;
k2=5.742728899;
t=0:0.001:10;   % time scalex 
initial_h=0; 
[t,h]=ode45( @rhs, t, initial_h); 
 
figure(1);
plot(t,h); 
xlabel('t'); ylabel('h'); 
figure(2);
plot(t,P1);
xlabel('t'); ylabel('P1');
 
    function dhdt=rhs(t,h) 
P1=k0+kM*(1+(ka-k1)/(k1-k2)*exp(-t/k1)+(ka-k2)/(k2-k1)*exp(-t/k2));
P2=P0+rho*g*h;
Qin=CV1*(P1-P2);
Qout=CV2*(P2-P3);
       dhdt = (Qout-Qin)/Area;
    end 
end
This is my code. When generating the plots of P1 versus time, it should look like a cubic function. However, it does not yield any results.
Also, the plot of h vs t is not accurate.
How do I fix this? Something seems off with the final result.
The ODE I am attempting to solve is
dh/dt=(Qout-Qin)/Area
Qout=CV2*(P2-P3);
Qin=CV1*(P1-P2);
For the Qout and Qin functions, P1 is dependent on time and P2 is dependent on h
P1=k0+kM*(1+(ka-k1)/(k1-k2)*exp(-t/k1)+(ka-k2)/(k2-k1)*exp(-t/k2));
P2=P0+rho*g*h;
For some reason matlab doesn't recognize this part of the code? I don't know how to fix this or if it needs to be moved somewhere else.

Best Answer

Since ‘rhs’ is a function nested inside ‘first_oder_ode’, it will inherit all the variables in the ‘first_oder_ode’ workspace. The reverse is not true, so you have to calculate ‘P1’ separately outside ‘rhs’ and after your ode45 call.
Your code with those changes: 
function first_oder_ode 
 
% SOLVE  dx/dt = -3 exp(-t). 
% initial conditions: x(0) = 0 
CV1=110*0.1336806;
CV2=110*0.1336806;
Area=15; % cross sectional area of tank, m^2
P0=1.0*10^5;
P3=1.1*10^5;
rho=1000;
g=9.8;
k0=10.09782398;
kM=423463.7642;
ka=8.869925491;
k1=0.633075579;
k2=5.742728899;
t=0:0.001:10;   % time scalex 
initial_h=0; 
[t,h]=ode45( @(t,h) rhs(t,h) , t, initial_h); 
P1=k0+kM*(1+(ka-k1)/(k1-k2)*exp(-t/k1)+(ka-k2)./(k2-k1).*exp(-t/k2));
figure(1);
plot(t,h); 
xlabel('t'); ylabel('h'); 
figure(2);
plot(t,P1);
xlabel('t'); ylabel('P1');
    function dhdt=rhs(t,h) 
P1=k0+kM*(1+(ka-k1)/(k1-k2)*exp(-t/k1)+(ka-k2)/(k2-k1)*exp(-t/k2));
P2=P0+rho*g*h;
Qin=CV1*(P1-P2);
Qout=CV2*(P2-P3);
       dhdt = (Qout-Qin)/Area;
    end 
end
Then just call it from a script or your Command Window as: 
first_oder_ode
to run it.
Related Question