How should I throttle my rocket to reach highest altitude? "Real world" problem.

Suppose we want to launch a rocket equipped with an engine which can be throttled as we prefer. Suppose also that the amount of fuel burnt per time is directly proportional to the upward force the engine exerts. How should we throttle the engine to reach the highest altitude, taking into account atmospheric drag?


Physical formulation.
Let $m(t) = m_d + m_f(t)$ be the total mass of the rocket, where $m_d$ is its dry mass and $m_f(t)$ the mass of the remaining fuel at time $t$. Our engine provides a force directly proprtional to the amount of fuel burnt; so the force exerted at the moment $t$ will be
$$ F = -\eta \dfrac{dm_f}{dt},$$
where $\eta$ is a positive constant which depends on the efficiency of the engine.
For the drag force we use the fluid-dynamics formulation:
$$ F_d = k(x) \dot x(t)^2, $$
where $x$ is our height and $k$ a coefficient which takes into account the change of density of the air with altitude.
So by Newton's law we have
$$\ddot x = \dfrac{F - F_d - mg}{m} = \dfrac{-\eta \dot{m_f} - k(x) \dot x ^2}{m_d + m_f}-g.$$
The problem has therefore become the following:
Find a non-negative function $\theta(t)$ (the throttle, i.e. $\theta = - \dot{m_f}$) such that
$$ \int_0^{\infty} \theta = m_f(0), \qquad \text{(we burn all the fuel, not more nor less)} $$
and such that the trajectory $x(t)$ solving the differential equation
$$ \ddot x(t) =  \dfrac{\eta \theta(t) - k(x(t)) \dot x(t)^2}{m_d + m_f(0) - \int_0^t \theta}-g$$
attains somewhere maximum value, that is for any other $\tilde \theta$ satsfying the previous requirements, its trajectory $\tilde x$ is such that $\sup_{t \in [0, \infty)} \tilde x(t) \leq \sup_{t \in [0, \infty)} x(t)$.

Mathematical formulation and recap.
Given $g, \eta, C>c \in (0, \infty), k : [0, \infty) \rightarrow [0,\infty) $ a decreasing function, let
\begin{align}
S_{c,C} := \Bigg\lbrace f \in L^1([0, \infty)) \text{ s.t. } & \int_0^{\infty} f = c ,\\
 & f \geq 0 \text{ a.e.} ,\\
& \exists x_f(t) \text{ s.t. } \ddot x_f(t) =  \dfrac{\eta f(t) - k(x_f(t)) \dot x_f(t)^2}{C - \int_0^t f}-g \Bigg\rbrace.
\end{align}
Find $\theta \in S_{c,C}$ s.t. $$\forall \tilde \theta \in S_C, \sup_{t \in [0, \infty)}  x_{\tilde \theta}(t) \leq \sup_{t \in [0, \infty)} x_{\theta}(t).$$
 A: This is an optimal control problem, so I will use the rules of optimal control.
First, we represent the state space equations. Also we take the total mass as a state and amount of fuel burnt as the input control. So we have:
\begin{cases} \tag{1}
\dot{x}_1=x_2 \\
\dot{x}_2=\frac{\eta \theta-k(x_1)x_2^2}{x_3}-g \\
\dot{x}_3=-\theta
\end{cases}
with these boundary conditions:
\begin{cases} \tag{2}
x_1(0)=0\\
x_2(0)=0\\
x_3(0)=m_d+m_f(0)\\
x_3(t_f)=m_d
\end{cases}
We want to find function $\theta(t)$ that minimize the following cost function:
$$J(\theta)=h(x(t_f),t_f)+\int_0^{t_f}g(x,\theta,t)dt=\int_{0}^{t_f}(-x_2)dt\tag{3}$$ 
In the above equations the terminal time $t_f$ is unknown and must be detemined.
For finding the solution, we define the Hamiltonian of this system:
\begin{eqnarray*}\tag{4}
H(t,\theta,x,p)&=&g+p^Ta\\
&=&-x_2+p_1x_2+p_2(\frac{\eta \theta-k(x_1)x_2^2}{x_3}-g)-p_3\theta
\end{eqnarray*}
where $a$ is the vector of states diffrential equations and $(p_1,p_2,p_3)$ are costate variables.
Using the calculus of variations, the necessary conditions for optimality will be find:
\begin{cases}\tag{5}
\dot{x}=\frac{\partial H}{\partial p}\\
\dot{p}=-\frac{\partial H}{\partial x}\\
0=\frac{\partial H}{\partial \theta}
\end{cases}
and we get this boundary condition:
\begin{eqnarray*}\tag{6}
&&\Biggl(\frac{\partial h}{\partial x}(x(t_f),t_f)-p(t_f) \Biggl)\delta x_f\\
&&+ \Biggl(H(x(t_f),p(t_f),\theta(t_f),t_f)+\frac{\partial h}{\partial t}(x(t_f),t_f) \Biggl)\delta t_f  \\
&& =0
\end{eqnarray*}
Using equations $(5)$ we can derive states and costates differential equations:
\begin{cases} \tag{7}
\dot{x}_1=x_2\\
\dot{x}_2=\frac{\eta \theta-k(x_1)x_2^2}{x_3}-g \\
\dot{x}_3=-\theta\\
\dot{p}_1=\frac{p_2x_2^2}{x_3}\frac{\partial k}{\partial x_1}\\
\dot{p}_2=1-p_1+2\frac{p_2k(x_1)x_2}{x_3}\\
\dot{p}_3=-\frac{p_2}{x_3^2}(\eta \theta-k(x_1)x_2^2)
\end{cases}
and from $0=\frac{\partial H}{\partial \theta}$ in $(5)$ we get:
$$\frac{p_2}{p_3}=\frac{x_3}{\eta} \tag{8}$$
Also from equation $(6)$ these final conditions will be find:
\begin{cases}\tag{9}
p_1(t_f)=0\\
p_2(t_f)=0
\end{cases}
and:
\begin{equation}\tag{10}
H(t_f)=0
\end{equation}
Now, we have seven unknown functions $\Bigl(x_1(t),x_2(t),x_3(t),p_1(t),p_2(t),p_3(t),\theta(t)\Bigl)$ with six differential equations in $(7)$ plus one constraint equation in $(8)$. For solving these equation we need sufficient boundary conditions. Using four boundary conditions in equation $(2)$ and two final conditions in equation $(9)$ all needed boundary condition will be known. 
We not used equation $(10)$ yet. This equation will be used for finding the terminal time $t_f$. 
Finally, for solving boundary value problem in $(7)$ and because of nonlinearity, we must use numerical methods such as shooting method.
A: I thought this question was interesting and I didn't want to do any proper work this afternoon so I made a simple model to find out what would happen. My matlab code is at the end of the question.
So far I've tested three cases and considered changing the initial thrust and adding a linear increase in the thrust for each case. The thrust is given as a fraction of the total fuel burnt per second.
Obviously the numbers I've used aren't incredibly realistic or sensible but they give an initial idea of what's going on. The color on the graphs represents the height obtained using that thrust.
So far I've only considered constant drag coefficients. Firstly I tried k=0 (no drag)

As expected this gives the greatest height when you burn fuel as fast as you can. This is mainly just a test to check I haven't done anything too wrong.
If we choose k=1 we get

Interestingly here we see the opposite effect. You want the lowest and longest thrust that will get you off the ground as the air resistance is the dominant form of drag.
Finally an intermediate case of k=0.1

Here there is a balance between air resistance and gravity and an intermediate value of thrust is favoured. 
Interestingly (at least to me), it was not advantageous to throttle up after launch. I suspect this may occur if variable k is considered.
Here is my matlab code. Should run in octave too.
function[x]=rocket_test(thrustparam)

g=9.81;

md=1000;
mf=1000;  
k=0.1;
nu=1000;

dt=0.001;

v=0;
x=0;
t=0;

c=thrustparam(1)*mf;
p_1=thrustparam(2)*mf;

count=1;

while v>=0                  %want to find max height keep going until you start descending.

    throttle=throttle_set(t, c, p_1);

    if mf-throttle*dt<0    %if you would burn more fuel than left burn all the fuel  left;
        throttle=mf/dt;
    end

m=md+mf;
Fd=k*v^2+m*g;
Ft=nu*throttle;

F=Ft-Fd;

a=F/m;

x=x+v*dt+0.5*a*dt^2;
v=v+a*dt;
mf=mf-throttle*dt;
t=t+dt;

%store(count, :)=[t, x, v, a, throttle, Fd, Ft, F, m];
%count=count+1;

end
end

function [throttle]=throttle_set(t, c, p_1)

throttle=c+p_1*t;

if throttle<0
    throttle=0;
end

end

