  We combine the parameterization method for invariant manifolds with the
finite element method for elliptic PDEs,to obtain a new computational framework
for high order approximation of invariant manifolds attached to unstable
equilibrium solutions of nonlinear parabolic PDEs. The parameterization method
provides an infinitesimal invariance equation for the invariant manifold, which
we solve via a power series ansatz. A power matching argument leads to a
recursive system of linear elliptic PDEs -- the so-called homological equations
-- whose solutions are the power series coefficients of the parameterization.
The homological equations are solved recursively to any desired order using
finite element approximation. The end result is a polynomial expansion for a
chart map of the manifold, with coefficients in an appropriate finite element
space. We implement the method for a variety of example problems having both
polynomial and non-polynomial nonlinearities, on non-convex two-dimensional
polygonal domains (not necessary simply connected), for equilibrium solutions
with Morse indices one and two. We implement a-posteriori error indicators
which provide numerical evidence in support of the claim that the manifolds are
computed accurately.
