In this paper we describe new B-spline Gaussian collocation software for solving two-dimensional
parabolic partial differential equations (PDEs) defined over a rectangular region. The numerical solution is represented as a bi-variate piecewise polynomial (using a tensor product
B-spline basis) with time-dependent unknown coefficients. These coefficients are determined by
imposing collocation conditions: the numerical solution is required to satisfy the PDE
and boundary conditions at images of the Gauss points mapped onto certain subregions of the spatial domain.
This leads to a large system of time-dependent
differential algebraic equations (DAEs) which is solved using the DAE solver, DASPK. We provide numerical results in which we use the new software, called BACOL2D, to solve three test problems.