  One goal in Bayesian machine learning is to encode prior knowledge into prior
distributions, to model data efficiently. We consider prior knowledge from
systems of linear partial differential equations together with their boundary
conditions. We construct multi-output Gaussian process priors with realizations
in the solution set of such systems, in particular only such solutions can be
represented by Gaussian process regression. The construction is fully
algorithmic via Gr\"obner bases and it does not employ any approximation. It
builds these priors combining two parametrizations via a pullback: the first
parametrizes the solutions for the system of differential equations and the
second parametrizes all functions adhering to the boundary conditions.
