  This paper introduces general methodologies for constructing closed-form
solutions to linear constant-coefficient partial differential equations (PDEs)
with polynomial right-hand sides in two and three spatial dimensions.
Polynomial solutions have recently regained significance in the development of
numerical techniques for evaluating volume integral operators and also have
potential applications in certain kinds of Trefftz finite element methods. The
equations covered in this work include the isotropic and anisotropic Poisson,
Helmholtz, Stokes, linearized Navier-Stokes, stationary advection-diffusion,
elastostatic equations, as well as the time-harmonic elastodynamic and Maxwell
equations. Several solutions to complex PDE systems are obtained by a potential
representation and rely on the Helmholtz or Poisson solvers. Some of the cases
addressed, namely Stokes flow, Maxwell's equations and linearized Navier-Stokes
equations, naturally incorporate divergence constraints on the solution. This
article provides a generic pattern whereby solutions are constructed by
leveraging solutions of the lowest-order part of the partial differential
operator (PDO). With the exception of anisotropic material tensors, no matrix
inversion or linear system solution is required to compute the solutions. This
work is accompanied by a freely-available Julia library,
\texttt{ElementaryPDESolutions.jl}, which implements the proposed methodology
in an efficient and user-friendly format.
