We propose efficient and parallel algorithms for the implementation of the high-order continuous time Galerkin method for dissipative and wave propagation problems. By using Legendre polynomials as shape functions, we obtain a special structure of the stiffness matrix which allows us to extend the diagonal Pad\'e approximation to solve ordinary differential equations with source terms. The unconditional stability, $hp$ error estimates, and $hp$ superconvergence at the nodes of the continuous time Galerkin method are proved. Numerical examples confirm our theoretical results.