The leap-frog algorithm is able to conserve to a certain extent the energy of a system, which flucutates as a cosine around a stable value. Is this true if we apply the algorithm to a n-body gravitational system? Does the energy of the system stay between two values or grows with time?
No, there are no guarantees.
The perturbation term for the energy is an expression in the derivatives of the Hamiltonian. For gravitational simulations, this Hamiltonian, more specifically the potential energy, has singular configurations, when the positions of two bodies coincide. This leads to non-reversible energy and momentum changes if two bodies get too close.
What "too close" is also depends on the step size. It should be easy to envision that a "fly-by" in 3 or so Verlet steps has not enough sampling points for the (locally) symmetrical nature of the force field to ensure the associated symmetry in the numerical solution.
This is not a problem in simulating a solar system with mostly circular orbits. However, in systems with bodies of an equal mass scale close encounters appear to be the norm rather than the exception, so any fixed-step method will eventually fail in keeping the energy constant.