MATLAB: Solving a coupled differential equations using ODE45

differential equationsMATLABode45runge kutta

I've a simple question about solving a differential equation using ode45. my only issue is that there is another differential equation in it and im not sure if what i am doing is wrong here is my code and im also attaching my differential functions.
.The initial condition for work done is w0=0.l i cant just write [th,W] = ode45(@work, range,w0)[gives me all 0s as output] because it is expressed in terms of function P(which also has its own initial condition P1) right ? i just need help in understanding how to write the function for work.
Thanks in advance.
P.S. I am aware that the pressure and work functions are similar. i just seperated them for troubleshooting sakes.
PFA the equations 
global  theta r bore stroke P1 T1 y  shaft rod epsilon V1 n a theta_s theta_d Qin  range
[th,P] = ode45(@pressure,range,P1);
plot(th,P)
function [V,dV] = vol(theta)
    global r ; 
    V = (1/r)  + ((r-1)*(1-cos(theta)))/2*r ;
    dV = ((r-1)/2*r )*sin(theta);
end
function [dxb] =burnfrac(theta)
    global theta_s theta_d n a ;
    xb = 1- exp(-a*(((theta-theta_s)*(1/theta_d))^n ));
    if theta>theta_s && theta<(theta_s+theta_d) 
        dxb= n*(a/theta_d)*(1-xb)*((theta-theta_s)/theta_d)^(n-1);
    else
        dxb=0;
    end
    
end
function [dP]= pressure(theta,P) 
    global Qin y ;
    [V ,dV] = vol(theta);
    dxb=burnfrac(theta);
    dQ= Qin *dxb ;
    dP= - y* P/V * dV + ((y-1)/V) * dQ;
    dW= P*dV;
   
   
end
function [dW]= work(theta,P) 
    global Qin y ;
    [V ,dV] = vol(theta)
    dxb=burnfrac(theta);
    dQ= Qin *dxb ;
    dP= - y* P/V * dV + ((y-1)/V) * dQ
    dW= P*dV
   
   
end

Best Answer

You have to put both ODEs together and integrate the pair at once. Let's look at an example equation:
That pair of equations can be simultaneosly integrated if we put them into one ode-function:
function dQdtdPdt = odePair(t,QP)
P = QP(2);
Q = QP(1); % for readability
dQdt = -P+f(t); % whatever you need for f...
dPdt = -P*Q + f;
dQdtdPdt = [dQdt;dPdt];
end
That function you now integrate with ode45 (for example):
Q0P0 = [12,32];
t_span = [0 34];
QP = ode45(@(t,QP) odePair(t,QP),t_span,Q0P0);
Which will give you both Q(t) and P(t).
HTH
Related Question