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.