  When numerical solution of elliptic and parabolic partial differential
equations is required to be highly accurate in space, the discrete problem
usually takes the form of large-scale and sparse linear systems. In this work,
as an alternative, for spatial discretization we provide a Matrix-Oriented
formulation of the classical Finite Element Method, called MO-FEM, of arbitrary
order $k\in\mathbb{N}$. On structured 2D domains (e.g. squares or rectangles)
the discrete problem is then reformulated as a Sylvester matrix equation, that
we solve by the reduced approach in the associated spectral space. On a quite
general class of domains, namely normal domains, and even on special surfaces,
the MO-FEM yields a multiterm Sylvester matrix equation where the additional
terms account for the geometric contribution of the domain shape. In
particular, we obtain a sequence of these equations after time discretization
of parabolic problems by the IMEX Euler method. We apply the matrix-oriented
form of the Preconditioned Conjugate Gradient (MO-PCG) method to solve each
multiterm Sylvester equation for MO-FEM of degree $k=1,\dots,4$ and for the
lumped $\mathbb{P}_1$ case. We choose a matrix-oriented preconditioner with a
single-term form that captures the spectral properties of the whole multiterm
Sylvester operator. For several numerical examples, we show a gain in
computational time and memory occupation wrt the classical vector approach
solving large sparse linear systems by a direct method or by the vector PCG
with same preconditioning. As an application, we show the advantages of the
MO-FEM-PCG to approximate Turing patterns with high spatial resolution in a
reaction-diffusion PDE system for battery modeling.
