I am trying to set up a code to calculate the equation of motion of a rocket using ode45

My equation for the EOM are

dv/dt= T/m - D/m - g*sin(gam): at=tangential acceleration, T=thrust, D=drag, g=acceleration of gravity, gam= path angle

dgam/dt= -1/v*(g-(v^2/(Re+h)))*cos(gam):  v=velocity, Re=radius of earth, h= altitude

dx/dt=Re/(Re+h)*v*cos(gam): x is downrange distance

dh/dt=v*sin(gam):

I have been watching videos, but I am not sure how to use ode45 with multiple functions, I would appreciate the help.

[Torsten]

Name v=y(1), gam=y(2), x=y(3) and h=y(4) and set up the ODE system for ODE45 like in the examples of the documentation.

[Jose Camorlinga Medina]

I was trying to follow that example, this is what I have so far

function dy = trajectory(t,y)

dy = zeros(4,1);

mf = mo-mpi; % Mass after propellent runs out (kg)

n = mo/mf; % Initial to final mass ratio ()

% time

t0 = 0; % Initial time (s)

tb = mf/md; % Burnout time (s)

tf = tb;

tspan = [t0,tf];

% Initial condition

v0 = 0; % Initial velocity (m/s)

gamma0 = 89.85*deg; % Initial flight path angle (rad)

x0 = 0; % Initial downrange distance (km)

h0 = 0; % Initial altitude (km)

vD0 = 0; % Initial value of velocity loss due to drag (m/s)

vG0 = 0; % Initial value of velocity loss due to gravity (m/s)

f0 = [v0; gamma0; x0; h0; vD0; vG0]; % Initial condition vector

% Equations of motion

dy(1) = (Ftv/m)-(D/m)-g*sin(y(2)); % Acceleration along path (m/s^2)

dy(2) = -(1/y(1))*(g-(y(1)^2/(Re+y(4))))*cos(y(2)); % Rate of change of path angle (rad)

dy(3) = Re/(Re+y(4))*y(1)*cos(y(2)); % Downrange distance (m)

dy(4) = y(1)*sin(y(2)); % Altitude (m)

[t,y] = ode45('trajectoty',tspan,f0)