syms w Ao q Ap Pa k L g t y(t) P(t) m(t) V c Y
Dy = diff(y);
D2y = diff(Dy);
DP = diff(P);
Dm = diff(m);
V = L+Dy;
Eq1 = D2y == -g+(P-Pa)*Ap
Eq2 = DP == k*P*((1/m) *(Dm) - (k/V)*(Dy))
Eq3 = Dm == m/V*(q-c*Ao*(Pa/P)^(1/k)*(2*k/(k-1) * (P*V/m)*(1-(Pa/P)^((k-1)/k)))^0.5)
[VF,Sbs] = odeToVectorField(Eq1, Eq2, Eq3)
Sample01 = matlabFunction(VF, 'Vars',{t,Y,k,q,L,Ao,c,Pa,Ap,g})
Sample01 = @(t,Y,k,q,L,Ao,c,Pa,Ap,g)[-(k.^2.*Y(1).*Y(3)-k.*q.*Y(1)+Ao.*c.*k.*Y(1).*(Pa./Y(1)).^(1.0./k).*sqrt((k.*(L+Y(3)).*(Pa-Y(1).*(Pa./Y(1)).^(1.0./k)).*(Pa./Y(1)).^(-1.0./k).*-2.0)./(Y(4).*(k-1.0))))./(L+Y(3));Y(3);-g+Ap.*Y(1)-Ap.*Pa;((q-Ao.*c.*(Pa./Y(1)).^(1.0./k).*sqrt((k.*(L+Y(3)).*(Pa-Y(1).*(Pa./Y(1)).^(1.0./k)).*(Pa./Y(1)).^(-1.0./k).*-2.0)./(Y(4).*(k-1.0)))).*Y(4))./(L+Y(3))];
Pa=101000;
d=.05;
Wp=7;
q=2;
Ao=0.01;
c=.1;
L=0.75;
k=1.4;
Te=300;
A=0.25*pi*d^2;
Ap = 0.42;
g = 9.81;
mass=(Pa*L*A)/(287*Te);
tic
W=200;
tf=10;
Fs = .1;
Ts = 1/Fs;
tspan =0:Fs:tf;
y0=[0,.0,101000,mass];
[t,y]=ode45(@(t,Y)Sample01(t,Y,k,q,L,Ao,c,Pa,Ap,g),tspan,y0);
subplot(2,2,1);plot(t,y(:,1),'k');
subplot(2,2,2);plot(t,y(:,2),'k');