  A newly developed weak Galerkin method is proposed to solve parabolic
equations. This method allows the usage of totally discontinuous functions in
approximation space and preserves the energy conservation law. Both continuous
and discontinuous time weak Galerkin finite element schemes are developed and
analyzed. Optimal order error estimates in both H^1 and L^2 norms are
established. Numerical tests are performed and reported.
