Evolution of a matter density under Newtonian gravity

Related evolution of current density and EM field.
This is my own question and answer in PSE.
QUESTION
I want to examine the physics of gravity. This subject has two parts: how the gravitational field influences the behavior of matter, and how matter determines the gravitational field. I want to understand first Newtonian gravity, where these two elements consist of the expression for the acceleration of a body in a gravitational potential $Φ$ , $$ \mathbf{a} = - \nabla \Phi, $$ and Poisson’s differential equation for the potential in terms of the matter density $ρ$ and Newton’s gravitational constant $G$ : $$ \nabla^2 \Phi = 4 \pi G \rho. $$ But the point is, the matter density is influenced itself by gravitation so, indeed we have $ρ = ρ (t, x)$ . Given $ρ (0, x)$ , how do I compute $ρ (t, x)$ ?

ANSWER
I have worked out my own answer. I would love someone tell me if this is all right or not.

In order to compute $ρ (t, x)$ given an initial density distribution $ρ (0, x)$ , we must take into account how the density evolves over time due to gravitational effects. This evolution can be described by attaching to Poisson equation above the continuity equation in fluid dynamics, combined with the Euler equation for the motion of the "infinitesimal particles".

More precisely:

Continuity Equation: With this we describe the conservation of mass and how the density $ρ$ changes over time:

\frac{\partial ρ}{\partial t} + \nabla \cdot (ρ v) = 0,

where $v (t, x)$ is the velocity field of the matter distribution. This equation ensures that mass is conserved as the matter moves.

Euler Equation (Momentum Conservation): This stands for the motion of the matter in response to forces, including the gravitational force:

\frac{\partial v}{\partial t} + (v \cdot \nabla) v = - \nabla Φ .

Here, $v (t, x)$ is the velocity field of the matter, and $Φ$ is the gravitational potential.

Solve the three PDEs simultaneously is probably impossible in most of the cases, so we proceed numerically as follows

Step 1: Use the initial density $ρ (0, x)$ and Poisson’s equation to compute the initial gravitational potential $Φ (0, x)$ .
Step 2: Solve the Euler equation to find the initial velocity field $v (0, x)$ of the matter.
Step 3: Using the velocity field and the continuity equation, update $ρ (t, x)$ over small time steps.
Step 4: Iterate this process. At each time step, update the gravitational potential $Φ (t, x)$ using Poisson’s equation, update the velocity field $v (t, x)$ using the Euler equation, and then update the density $ρ (t, x)$ using the continuity equation.