How to minimise jettisoned mass spent for a body with continuous but variable thrust?

The usual approach is to use optimal control. I will change the notations slightly. Let $x$ be the position, $v$ the velocity and $m$ the mass. I will also set the constant the exhaust relative speed $|v_e|=1$. You can combine the direction of the rocket and mass exhaust into a single real control parameter $u = -\frac{\dot m}m v_e$. The equations of motion are: $$ \dot x = v\quad \dot v = u-g $$ You want to minimise the total mass expelled (setting the initial mass to $1$ with no loss of generality) during the time interval $t\in(0,1)$, you will need to minimise the action: $$ S = \int_0^1|u|dt $$ under the constraint that $x(0) = 0,v(0) = v_0,x(1) = x_1$ and along the trajectory. You therefore identify the Lagrangian and compute the corresponding Hamiltonian: $$ L(x,v,u) = |u|\quad H(x,v,p_x,p_v) = \sup_up_xv+p_v(u-g)-|u| = \begin{cases} -\infty & p_v<-1\\ p_xv-p_vg & -1<p_v<1\\ +\infty & 1<p_v\end{cases} $$ Essentially, this is telling you that for optimal trajectories, $u=0$ or $u=\pm\infty$ so that you instantaneously eject a finite amount of mass. This makes sense, in terms of distance travelled, ejection $udt$ at time $t_*\in(0,1)$ will result distance $udt(1-t_*)$, so to make the most of it, you just need to eject everything at the beginning $t=0$ instantaneously. This amounts to effectively tuning the initial velocity just to reach the final position in the desired time under free fall. If you also want to additionally control the final velocity (consider landing for example), then you would additionally need to add a final instantaneous boost accordingly.

For practical purposes, you therefore need to refine your objective function. Everything will depend on how you penalise high mass ejection rate $|u| \to +\infty$. A simple tweak is to impose a maximum value $|u|\leq U$. Another analytically motivated amendment is to replace it by add a quadratic cost, i.e. a soft boundary, $u^2/2$. For the former choice, the new Hamiltonian is: $$ H(x,v,p_x,p_v) = \sup_{u\in[-U,U]}p_xv+p_v(u-g)-|u| = \begin{cases} p_xv-p_v+U(|p_v|-1) & p_v\not\in(-1,1)\\ p_xv-p_vg & p_v\in(-1,1)\end{cases} $$

Applying Pontryaguin's principle, you get the equations of motion: $$ \dot x = v\quad \dot v = \begin{cases} -g & p_v\in(-1,1)\\ U\operatorname{sgn}(p_v)-g & p_v\not\in(-1,1)\end{cases}\\ \dot p_x = 0\quad \dot p_v = -p_x $$ The trajectories therefore comprise of on/off states where the rocket is maximally thrusting or not at all respectively. By tuning the piecewise parabolas, you can match the beginning and final desired state. Note that you now have an issue of attainability if $U$ is not large enough to accommodate for the desired distance and velocity.

Hope this helps.