  We develop a formalism for General Relativistic N-body simulations in the
weak field regime, suitable for cosmological applications. The problem is kept
tractable by retaining the metric perturbations to first order, the first
derivatives to second order and second derivatives to all orders, thus taking
into account the most important nonlinear effects of Einstein gravity. It is
also expected that any significant "backreaction" should appear at this order.
We show that the simulation scheme is feasible in practice by implementing it
for a plane-symmetric situation and running two test cases, one with only cold
dark matter, and one which also includes a cosmological constant. For these
plane-symmetric situations, the deviations from the usual Newtonian N-body
simulations remain small and, apart from a non-trivial correction to the
background, can be accurately estimated within the Newtonian framework. The
correction to the background scale factor, which is a genuine "backreaction"
effect, can be robustly obtained with our algorithm. Our numerical approach is
also naturally suited for the inclusion of extra relativistic fields and thus
for dark energy or modified gravity simulations.
