The technique to do this is similar to that used in constructing stellar models. You know some of the properties of your object - in this case, it look like the mass and radius. You want to figure out the internal structure of the planet, including the central density and pressure, as well as the density and pressure profiles as functions of radius. The best way to do this is as follows:
- Determine a relationship between density and pressure - what we call an equation of state. We write this as $\rho = f(P)$.
- Make a guess as to the central pressure - and therefore the central density, using the above relation.
- Numerically integrate the differential equations determining the structure of the planet. Start in the core, and go outwards.
- Once your equations indicate that you've reached the surface of the planet, see what the predicted mass and radius are.
- Given the result of step 4, adjust your guesses for the central pressure and density, and repeat steps 3 and 4.
The two equations you need, besides the equation of state, are the equation of hydrostatic equilibrium $(1)$ and the equation of mass continuity $(2)$:
$$\frac{dP}{dr} = -\frac{GM_r\rho}{r^2}\tag{1}$$
$$\frac{dM}{dr} = 4\pi r^2\rho\tag{2}$$
Here, $P$ is pressure, $\rho$ is density, $r$ is radius, $G$ is the gravitational constant, and $M_r$ is the mass contained within a radius $r$:
$$M_r = \int_0^r 4\pi r'^2\rho(r')dr'$$
Let's say that at a given radius $r_i$, you know the pressure, density, and mass encapsulated within that radius - $P_i$, $\rho_i$, and $M_{r,i}$. Assume that your integration uses some step size $\Delta r$. Then the values at the radius $r_{i+1}$ are
$$r_{i+1} = r_i + \Delta r$$
$$P_{i+1} = P_i + \Delta P = p_i + \frac{dP}{dr}\Delta r = P_i - \frac{GM_{r,i}\rho_i}{r_i^2}\Delta r$$
$$\rho_{i+1} = f(P_{i+1})$$
$$M_{r,i+1} = M_i + \Delta M = M_i + \frac{dM}{dr}\Delta r = M_i + 4\pi r_i^2\rho_i\Delta r$$
Continue doing this until you reach the surface of the planet, i.e. at $P = 0$, then repeat steps 3 to 4 with a new guess until you're satisfied with the model's predicted mass and radius.
I implemented this procedure in my answer to another question, although I only had to do steps 3 to 4 once - I had a reasonable central pressure all figured out. I used an equation of state from Seager et al. 2008, designed for rocky planets.