Yes, you can derive the equations of motion for
$N$ particles (with masses $m_i$ at positions $\mathbf{x}_i(t)$, for $i=1, ... N$)
and for the gravitational field ($\Phi(\mathbf{x},t)$, for $\mathbf{x}\in\mathbb{R}^3$)
from a common Lagrangian.
However, we need to carefully distinguish between
the discrete positions $\mathbf{x}_i$ of the particles
and the continuous position $\mathbf{x}$ of the field.

The complete Lagrangian looks like this:

$$L[\mathbf{x}_i,\dot{\mathbf{x}}_i,\Phi,\nabla\Phi]=
  \sum_{i=1}^N \left(\frac{1}{2}m_i\dot{\mathbf{x}}_i(t)^2
  -m_i \Phi(\mathbf{x}_i(t))\right)
  -\frac{1}{8\pi G}\int d^3 x\ (\nabla\Phi(\mathbf{x},t))^2  \tag{1}$$
It has 3 components:
* a kinetic energy for each particle,
* an interaction energy between each particle and the gravitational field,
* the energy of the gravitational field.

We can define the density field made up by the particles as
$$\rho(\mathbf{x},t)=\sum_{i=1}^N m_i\ \delta^3(\mathbf{x}-\mathbf{x}_i(t))  \tag{2}$$
which has Dirac-like peaks where the particles are.

Using this density we can rewrite the Lagrangian (1). It is equivalent to:

$$L[\mathbf{x}_i,\dot{\mathbf{x}}_i,\Phi,\nabla\Phi]=
  \sum_{i=1}^N \frac{1}{2}m_i\dot{\mathbf{x}}_i(t)^2
  +\int d^3 x\underbrace{\left(
  -\rho(\mathbf{x},t)\Phi(\mathbf{x},t)-\frac{1}{8\pi G}(\nabla\Phi(\mathbf{x},t))^2
  \right)}_{\mathcal{L[\Phi,\nabla\Phi]}}  \tag{3}$$

From Lagrangian (1) or (3) we get the equations of motion.

From Lagrangian (1) the Euler-Lagrange equations
with respect to $\mathbf{x}_i$ give Newton's law:
$$\frac{d}{dt}\left(\frac{\partial L}{\partial\dot{\mathbf{x}}_i}\right)
  =\frac{\partial L}{\partial\mathbf{x}_i}$$
$$m_i\ddot{\mathbf{x}}_i=-m_i\nabla\Phi(\mathbf{x}_i)  \tag{4}$$

And from Lagrangian (3) the Euler-Lagrange equation
with respect to $\Phi(\mathbf{x},t)$
gives Poisson's equation:
$$\nabla\left(\frac{\partial{\mathcal{L}}}{\partial(\nabla\Phi)}\right)
  =\frac{\partial\mathcal{L}}{\partial\Phi}$$
$$-\frac{1}{4\pi G}\nabla^2\Phi=-\rho  \tag{5}$$