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.