This project is concerned with the development of implementing new matrix decomposition algorithms for solving the systems of linear algebraic equations arising when orthogonal spline collocation (that is, spline collocation at Gauss points) is applied to separable, second order, linear problems on rectangular regions. These algorithms have application in the solution of nonseparable boundary value problems, problems on general regions, and time dependent problems. The success of matrix decomposition algorithms depends on knowledge of the eigenvalues and eigenvectors of matrices arising in corresponding two-point boundary value problems in one space dimension. Knowledge of such an eigensystem with respect to one coordinate variable reduces the original two-dimensional discrete problem to solving a collection of independent discrete two-point boundary value problems with respect to the other coordinate variable. A software package is developed implementing orthogonal spline collocation algorithms with polynomials of arbitrary order on nonuniform partitions in each coordinate direction. The package includes, as a special case, a fast transform solver for piecewise Hermite bicubic orthogonal spline collocation for Poisson's equation in various coordinate systems.