clc;
clear all;
global m0 g0 T A Cd rh0 H0 Re hgr_turn tf md
% Launch Site: Guiana Space Center
Alt = 1;              %[m] Alt above sea level

% VEGA Rocket
m_stage_gross = [95796, 25751,10948];% 1st, 2nd,3d
% First stage(Solid Fuel)
m_prop = 88365;       % [kg] Propellant mass
Isp    = 280 ;        % [s]  Specific impulse
d      = 3;           % [m]  Diameter
g0     = 9.81;        % [m/s^2] Constant at its sea-level value
m0  = 137000;         % [kg] Initial mass
A   = pi*d^2/4;       % [m^2]Frontal area
Cd  = 0.5 ;             % Drag coefficient,assumed to have the constant value
rh0 = 1.225;          % [kg/m^3]
H0 = 7500;            % [m] Density scale height
Re = 6378e3;          % [m] Earth's radius
hgr_turn = 200;       % [m] Rocket starts the gravity turn when h = hgr_turn
tburn = 106.8;        % [s] Fuell burn time, first stage
md = (m_prop)/tburn;  % [kg/s]Propellant mass flow rate
T   = md*(Isp*g0);    % [N] Thrust (mean)
mf = m0 - m_prop;     % [kg] Final mass of the rocket(first stage is empty)
t0 = 0;               % Rocket launch time
tf = t0 + tburn;      % The time when propellant is completely burned
%and the thrust goes to zero
t_range     = [t0,tf];  % Integration interval

% Launch initial conditions:
gamma0 = 89.5/180*pi;       % Initial flight path angle
v0 = 0;   % Velocity (m/s)  % Earth's Rotation considered in eq of motion.
x0 = 0;   % Downrange distance [km]
h0 = Alt; % Launch site altitude [km]
vD0 = 0;  % Loss due to drag (Velocity)[m/s]
vG0 = 0;  % Loss due to gravity (Velocity)[m/s]
state0   = [v0, gamma0, x0, h0, vD0, vG0];
% Solve initial value problem for ordinary differential equations
[t,state] = ode45(@RocketDynEq,t_range,state0) ;
v     = state(:,1)/1000;      % Velocity [km/s]
gamma = state(:,2)*180/pi;    % Flight path angle  [deg]
x     = state(:,3)/1000;      % Downrange distance [km]
h     = state(:,4)/1000;      % Altitude[km]
vD    = -state(:,5)/1000;     % Loss due to drag (Velocity)[m/s]
vG    = -state(:,6)/1000;     % Loss due to gravity (Velocity)[m/s]
plot(t,h,'b');
hold on;
grid on;
plot(t,h,'.b');
title('Rocket Dynamics,VEGA Rocket sub-orbital ascent profile(1st stage P80)');
xlabel('time[s]');
ylabel('Altitude[km]');
 text(80,5,'smallsats.org','Color',[0 0 1], 'VerticalAlignment','middle',...
	'HorizontalAlignment','left','FontSize',14 );

% VEGA Rocket: First Stage P80
fprintf('\n VEGA Rocket: First Stage P80\n')
fprintf('\n Propellant mass           = %4.2f [kg]',m_prop)
fprintf('\n Gross mass                = %4.2f [kg]',m_stage_gross(1))
fprintf('\n Isp                       = %4.2f [s]',Isp)
fprintf('\n Thrust(mean)              = %4.2f [kN]',T/1000)
fprintf('\n Initial flight path angle = %4.2f [deg]',gamma0*180/pi)
fprintf('\n Final speed               = %4.2f [km/s]',v(end))
fprintf('\n Final flight path angle   = %4.2f [deg]',gamma(end))
fprintf('\n Altitude                  = %4.2f [km]',h(end))
fprintf('\n Downrange distance        = %4.2f [km]',x(end))
fprintf('\n Drag loss                 = %4.2f [km/s]',vD(end))
fprintf('\n Gravity loss              = %4.2f [km/s]',vG(end))
fprintf('\n');

 VEGA Rocket: First Stage P80

 Propellant mass           = 88365.00 [kg]
 Gross mass                = 95796.00 [kg]
 Isp                       = 280.00 [s]
 Thrust(mean)              = 2272.67 [kN]
 Initial flight path angle = 89.50 [deg]
 Final speed               = 1.76 [km/s]
 Final flight path angle   = 80.11 [deg]
 Altitude                  = 71.51 [km]
 Downrange distance        = 54.55 [km]
 Drag loss                 = 0.05 [km/s]
 Gravity loss              = 1.04 [km/s]

VEGA Rocket sub-orbital ascent profile(1st stage P8

function dfdt = RocketDynEq(t,y)
global m0 g0 T A Cd rh0 H0 Re hgr_turn md
v  =  y(1);     % Velocity
gm =  y(2);     % Flight path angle
x  =  y(3);     % Downrange distance
h  =  y(4);     % Altitude
vD =  y(5);     % Velocity loss due to drag
vG =  y(6);     % Velocity loss due to gravity
% Equations of motion of a gravity turn trajectory
      m = m0 - md*t;  % Vehicle mass
% else
%     m = mf;          % Burnout mass
%     T = 0;           % No more thrust is generated
% end
g  = g0/(1 + h/Re)^2;          % Gravitational variation with altitude
rh = rh0*exp(-h/H0);            % Atmospheric density exponential model
D = 1/2 * rh*v^2 * A * Cd;      % Drag force

% Rocket starts the gravity turn when h = hgr_turn
if h <= hgr_turn % Vertical flight
    dv_dt  = T/m - D/m - g;
    dgm_dt = 0;
    dx_dt  = 0;
    dh_dt  = v;
    dvG_dt = -g;
else                                  % Gravity turn
    dv_dt  = T/m - D/m - g*sin(gm);
    dgm_dt = -1/v*(g - v^2/(Re + h))*cos(gm);
    dx_dt  = Re/(Re + h)*v*cos(gm) + 463*sin(gm);
    dh_dt  = v*sin(gm) + 463*cos(gm); % Adding earth's rotation speed
    dvG_dt = -g*sin(gm);              % Gravity loss rate [m/s^2]
end
    dvD_dt = -D/m;           % Drag loss rate [m/s^2]
dfdt = [ dv_dt,dgm_dt, dx_dt,dh_dt, dvD_dt, dvG_dt]';
return

Image credit:  http://www.arianespace.com/launch-services-vega/VEGAUsersManual.pdf