Analytic and numerical integration of the rocket equation yield different results [closed]

Background

The mass of a rocket as a function of time is: $$m = m_0 - \dot m t$$ Where $$m_0$$ is the initial mass of the rocket and $$\dot m$$ is the mass flow rate. At a certain time all the fuel has been used, so from that time onwards $$m = m_f$$, the mass of the rocket without fuel.

The speed is then defined by the following equation, where $$V_e$$ is the exhaust velocity and $$g$$ is constant gravity. $$v = V_e log({m_o \over m_o - \dot mt}) - gt$$

To obtain position let $$\dot m$$ be constant, the equation can be written as

$$v = V_e log({m_o \over m_o - kt}) - gt$$ Integration of that expression lets to the following, which corresponds with the expression contained in Orbital Mechanics by Curtis.

$$h = \frac{V_e}{\dot m}[(m_0 - \dot m t)log({m_o - \dot mt\over m_o })+ \dot mt] - \frac{1}{2}g t^2$$

Question

However, the plot of $$h$$ (Red) is different for what I get for evaluating $$v$$ and then doing a numerical integration (Blue). Both plots are similar until the time where the fuel runs out, then the difference is easily visible.

The analytic plot sounds more reasonable to me since a rocket needs multiple stages to achieve a high altitude, but the numerical method is harder to get wrong and is necessary to consider other factors such as drag. Can someone point me in the right direction? Am I doing something wrong?

Here is the code I used in MATLAB. The values used are an approximation to the Saturn V, just considering one stage.

t = 0:1:1600;

%Constants
mass_inital = 2970.000;
mass_final = 731.000;
flow = 13.120;
exhaust = 2.500;
g = 9.81/1000;

mass = mass_inital - flow .* t;
mass(mass< mass_final) = mass_final;

v = exhaust .* log(mass_inital./mass) - (g .* t);

c = cumtrapz(v);
h = (exhaust ./ flow).* (mass.*log(mass./ mass_inital)+ flow .* t) - (0.5 .*g .*(t.*t));

%Plots
plot(t,c);
hold on
plot(t,h);
ylim([0 max(c)]); closed as off-topic by Kyle Kanos, user191954, John Rennie, Jon Custer, ZeroTheHeroOct 17 '18 at 0:15

• This question does not appear to be about physics within the scope defined in the help center.
If this question can be reworded to fit the rules in the help center, please edit the question.

• I'm voting to close this question as off-topic because it's about debugging code and not physics. – Kyle Kanos Oct 15 '18 at 21:36
• "the numerical method is harder to get wrong" - that is one of the funnier things I've read in a while! People get numerical methods wrong all the time... – Jon Custer Oct 16 '18 at 19:19

Most likely the last term in the square brackets for h ($$\dot{m}t$$) is incorrect when fuel is exhausted and is imposing extra height.