  We describe a numerical approach to address the BKL conjecture that the
generic cosmological singularity is locally Mixmaster-like. We consider
application of a symplectic PDE solver to three models of increasing
complexity--spatially homogeneous (vacuum) Mixmaster cosmologies where we
compare the symplectic ODE solver to a Runge-Kutta one, the (plane symmetric,
vacuum) Gowdy universe on $T^3 \times R$ whose dynamical degrees of freedom
satisfy nonlinearly coupled PDE's in one spatial dimension and time, and U(1)
symmetric, vacuum cosmologies on $T^3 \times R$ which are the simplest
spatially inhomogeneous universes in which local Mixmaster dynamics is allowed.
