For a parabolic problem with maximal $L^p$-regularity, we prove that the time discretization by a linear multistep method or Runge-Kutta method preserves the maximal $L^p$-regularity uniformly in the stepsize if the method is A-stable (and satisfies minor additional conditions). In particular, the implicit Euler method, the Crank-Nicolson method, the second-order backward difference formula (BDF), and the Radau IIA and Gauss Runge-Kutta methods of all orders preserve maximal regularity. In general polygons and polyhedra, possibly nonconvex, the analyticity of the finite element heat semigroup in the $L^q$ norm, $1 \le q \le \infty$, and the maximal $L^p$-regularity of finite element solutions of parabolic equations are proved.