I want to see how much fuel it takes to send a rocket into space. By changing the initial fuel, after I pass 2150 tones, the speed vy becomes negative. Even by keeping the fuel under 2150 tones the rocket doesn't manage to fly so far. Is my implementation correct?The equations I am using are:\begin{align}  &\frac{d v_x}{dt} = \frac{\alpha \Theta_t u_x}{m(t)} - \frac{G M x}{(x^2 + y^2)^{3/2}}\\  &\frac{d v_y}{dt} = \frac{\alpha \Theta_t u_y}{m(t)} - \frac{G M y}{(x^2 + y^2)^{3/2}}\\  &\frac{dm}{dt} = -\alpha \Theta_t\end{align}which are from a book (link here). I have modified it slightly as follows\begin{align}  &a_x = \frac{F_{thrust_x} - m g_{gravity_x} - F_{drag_x}}{m}\\  &F_{gravity_x} = G \frac{Mx}{(x^2+y^2)^\frac{3}{2}}\end{align}Thrust is a constant until fuel runs out, then it is zero.The above is also valid for the y component by changing the subscripts and the $x$ to $y$ in the gravity function.Function:function [t,sol] = GetAcceleration_7 ( t,h,rocket,fuel,ic,m,ft)    % ic = [x,vx,y,vy];    sol = rungekuttastep_xu(t,h, fuel, ic, rocket,m,ft);endFunction:function sol = rungekuttastep_xu(t,h, fuel, w, R,m,F)    k1 = f(t,w,fuel, R,m,F);    k2 = f(t+h/2,w+h*k1/2,fuel, R,m,F);    k3 = f(t+h/2,w+h*k2/2,fuel, R,m,F);    k4 = f(t+h, w+h*k3, fuel, R,m,F);    sol = w+(k1 + 2*k2 + 2*k3 + k4)*h/6;endfunction w = f(t,w,fuel, R,m,F)    G = 6.673e-11; M = 5.98e24;    gx = (G*M*w(1))/(w(1)^2 + w(3)^2)^(3/2);    gy = (G*M*w(3))/(w(1)^2 + w(3)^2)^(3/2);    dx = 0.8 * w(2)^2; dy = 0.8 * w(4)^2;    angle = 55;    theta = 2*angle*pi/360;    Fx = F*cos(theta);    Fy = F*sin(theta);    w(1) = w(2);    if w(2) >= 0        w(2) = (Fx - m * gx - dx)/m;         else        w(2) = (Fx - m * gx + dx)/m;         end    w(3) = w(4);    if w(4) >= 0        w(4) = (Fy - m * gy - dy)/m;        else        w(4) = (Fy - m * gy + dy)/m;         endendMain Program:% Apollo% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Apollo.initialRocketMass = 750e3;                           % rocket massApollo.initialFuelMass = 150e3; %2150e3;                     % fuel massApollo.totalMass = Apollo.initialRocketMass...                    + Apollo.initialFuelMass;               % total massApollo.burnRate = 15e3;                                     % burn rateApollo.initialThrust = 34e6;                                % initial thrust%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%dt = 0.1;               % time steptstart = 0;             % start timetstop = 1000;           % upper limit on how long to simulateA = []; V = []; H = []; T = []; M = []; F = []; % empty vectorsk = 1;y0 = 6.5e6; vx = 0*dt; vy = 0;      % initial velocityx = 0; y = y0;          % initial heightM(1) = Apollo.totalMass;    % enter the loop with total mass    F(1) = Apollo.initialThrust;Y(1) = y;                   % enter the loop with initial positionm = M(1);ft = F(1);w = [x,vx,y,vy];for t = tstart:dt:tstop    if sqrt(w(1)^2+w(3)^2) >= y0 %ground is y0        % track variable mass        if ( t <= 0 )            m = Apollo.totalMass;            display('1')        elseif ( t > 0 && fuel == true )            m = Apollo.totalMass - Apollo.burnRate*t;            display('2')        else            m = Apollo.initialRocketMass;            display('3')        end        % track fuel consumption        if m > Apollo.initialRocketMass            fuel = true;            display('4')        else            fuel = false;            display('5')        end        % track thrust        if ( t >= 0 && Apollo.initialFuelMass >= Apollo.burnRate && fuel == true)            ft = Apollo.initialThrust;            display('6')        else            ft = 0;            display('7')        end        [t, w] = GetAcceleration_7( t,dt,Apollo,fuel, w,m,ft); % get current acceleration        v = sqrt((w(2))^2+(w(4))^2); % save speed        % telemetry for plotting        k = k + 1;          V(k)= v ; VX(k) = w(2); VY(k) = w(4);        X(k)= w(1) ; Y(k)= w(3) ;        M(k) = m; F(k) = ft;         T (k)= t ;         D(k) = sqrt(w(1)^2+w(3)^2);   % track max distance traveled    endend
Two stage rocket launch using Runge-Kutta 4 gives weird results
matlab;numerical analysis;runge kutta