The non-homogenous 1D Euler equations with self-gravity are given by

\frac{\partial }{\partial t} \begin{pmatrix} \rho \\ \rho u \\ E\end{pmatrix} + \frac{\partial }{\partial x} \begin{pmatrix} \rho u \\ \rho u^2 + p \\ u(E+p) \end{pmatrix} = \begin{pmatrix} 0 \\ - \rho \triangledown \phi_g \\ -\rho u \cdot \triangledown \phi_g \end{pmatrix}

which in short hand notation is

\vec{u}_t + \vec{F}(\vec{u})_x = \vec{S}(\vec{u})

where \vec{u} is vector of fluid variables, \vec{F} is their fluxes, and \vec{S} is the source-term vector.

To solve this system of equations, we can employ operator splitting, which is analagous to the procedure for splitting the higher dimensional Euler equations (as I did for 2D — see wiki page).

It is basically as follows, first solve the homogenous equations,

\vec{u}_t + \vec{F}(\vec{u})_x = 0

with initial condition for the grid,

\vec{u}(x,0) = \vec{u}_i ^0

over a time-step dt (found by usual CFL condition for upwind Godunov scheme). This gives the solution,

\vec{u}^{n+1/2}

Next, solve the equation for the source-term,

\frac{d}{dt}\vec{u} = \vec{S}(\vec{u})

with initial condition

\vec{u}(x,0) = \vec{u}^{n+1/2}

This gives the solution for the complete time-step,

\vec{u}^{n+1}

Given the ODE's for the 'source' step only involve equations for momentum and energy, we conclude that only u and E change over this step (density does not).

Thus, a schematic is as follows,

Hydro step \rightarrow \vec{u} = <\rho^{n+1/2}, u^{n+1/2}, E^{n+1/2}>~\rightarrow

Source step \rightarrow \vec{u} = <\rho^{n+1/2}=\rho^{n+1}, u^{n+1}, E^{n+1}>