Need Help Solving System of ODEs

close all; clear all; clc

I wanted to solve these 4 ode(s) and was trying to use ode45, but the output is not what I was expecting, the v is supposed to decrease but it was increasing, the gamma was also supposed to become more negetive but it remained constant. Any idea about where I must be going wrong? Thanks in advance!

g=9.81;
z0=10000;
R=6378000;
h0=(R*z0)/(R+z0);
for z=1:10000
    h(z)=(R.*z)./(R+z);
    rho0=0.4135;
    H=8500;
    rho=rho0.*exp((h0-h)./H);
end

v0=450;
gamma0=-20;
m0=2500;
T=15000;
Isp=300;

dodt=@(t,W)[(-g*sin(W(2)))-(T/W(4));(-g/W(1))*cos(W(2));W(1)*sin(W(2));-T/(g*Isp)];
[t,W]=ode45(dodt,[0,21.93],[v0,gamma0,z0,m0])

v=W(:,1);
gamma=W(:,2);
z=W(:,3);
m=W(:,4);

Answer 1

[Alan Stevens]

You haven't included the first term on the right-hand side of your equation for dV/dt, nor for dgamma/dt.

[Akrant Rai]

Thank you for the response.I forgot to mention that those terms are essentially 0 for the case I am working on and thus have been neglected. Also phi is 0 degrees everywhere and thus cos phi and sin phi have been treated as 1 and 0.

[Alan Stevens]

Try

dodt=@(t,W)[(-g*sin(W(2)))-(T/W(4));(-g/W(1))*cos(W(2))+T/(g*W(1));W(1)*sin(W(2));-T/(g*Isp)];

There was a term missing : T/(g*W(1)) , in the equation for dgamma/dt.

[Akrant Rai]

Did you mean T/(W(4)*W(1))? That term is 0 since phi and thus sin phi is 0. Thanks anyway though!

[Alan Stevens]

I did mean that - I misread the term! So, I don't see anything else wrong with the equations! You will need to check they are derived correctly, and that you have used the right data.

[Akrant Rai]

Alright. I will check the physics again to make sure that is right. Thank you for all the help!