  The energy balance of a partially molten rocky planet can be expressed as a
non-linear diffusion equation using mixing length theory to quantify heat
transport by both convection and mixing of the melt and solid phases. In this
formulation the effective or eddy diffusivity depends on the entropy gradient,
$\partial S/\partial r$, as well as entropy. First we present a simplified
model with semi-analytical solutions, highlighting the large dynamic range of
$\partial S/\partial r$, around 12 orders of magnitude, for physically-relevant
parameters. It also elucidates the thermal structure of a magma ocean during
the earliest stage of crystal formation. This motivates the development of a
simple, stable numerical scheme able to capture the large dynamic range of
$\partial S/\partial r$ and provide a flexible and robust method for
time-integrating the energy equation.
  We then consider a full model including energy fluxes associated with
convection, mixing, gravitational separation, and conduction that all depend on
the thermophysical properties of the melt and solid phases. This model is
discretised and evolved by applying the finite volume method (FVM), allowing
for extended precision calculations and using $\partial S/\partial r$ as the
solution variable. The FVM is well-suited to this problem since it is naturally
energy conserving, flexible, and intuitive to incorporate arbitrary non-linear
fluxes that rely on lookup data. Special attention is given to the numerically
challenging scenario in which crystals first form in the centre of a magma
ocean.
  Our computational framework is immediately applicable to modelling high melt
fraction phenomena in Earth and planetary science research. Furthermore, it
provides a template for solving similar non-linear diffusion equations arising
in other disciplines, particularly for non-linear functional forms of the
diffusion coefficient.
