We present a new time discretization scheme based on the continuous Galerkin Petrov
method of polynomial order 3 (cGP(3)-method) which is combined with a reduced numerical
time integration (3-point Gauß-Lobatto formula). The solution of the new approach can be
computed from the solution of the lower order cGP(2)-method, which requires to solve a
coupled 2 × 2 block system on each time interval, followed by a simple post-processing step,
such that we get the higher accuracy of 4th order in time in the standard L2-norm with
nearly the cost of the cGP(2)-method. Moreover, the difference of both solutions can be
used as an indicator for the approximation error in time. For the approximation in space we
use the nonparametric eQ3-element which belongs to a family of recently derived higher order
nonconforming finite element spaces and leads to an approximation error in space of order
4, too, in the L2-norm. The expected optimal accuracy of the full discretization error in the
L2-norm of 4th order in space and time is confirmed by several numerical tests. We discuss
implementation aspects of the time discretization as well as efficient multigrid methods for
solving the resulting block systems which lead to convergence rates being almost independent
of the mesh size and the time step. In our numerical experiments we compare different
higher order spatial and temporal discretization approaches with respect to accuracy and
computational cost.