  We propose a Bernoulli-barycentric rational matrix collocation method for
two-dimensional evolutionary partial differential equations (PDEs) with
variable coefficients that combines Bernoulli polynomials with barycentric
rational interpolations in time and space, respectively. The theoretical
accuracy $O\left((2\pi)^{-N}+h_x^{d_x-1}+h_y^{d_y-1}\right)$ of our numerical
scheme is proven, where $N$ is the number of basis functions in time, $h_x$ and
$h_y$ are the grid sizes in the $x$, $y$-directions, respectively, and $0\leq
d_x\leq \frac{b-a}{h_x},~0\leq d_y\leq\frac{d-c}{h_y}$. For the efficient
solution of the relevant linear system arising from the discretizations, we
introduce a class of dimension expanded preconditioners that take the advantage
of structural properties of the coefficient matrices, and we present a
theoretical analysis of eigenvalue distributions of the preconditioned
matrices. The effectiveness of our proposed method and preconditioners are
studied for solving some real-world examples represented by the heat conduction
equation, the advection-diffusion equation, the wave equation and telegraph
equations.
