  We propose a collocation method based on multivariate polynomial splines over
triangulation or tetrahedralization for the numerical solution of partial
differential equations. We start with a detailed explanation of the method for
the Poisson equation and then extend the study to the second-order elliptic PDE
in non-divergence form. We shall show that the numerical solution can
approximate the exact PDE solution very well. Then we present a large amount of
numerical experimental results to demonstrate the performance of the method
over the 2D and 3D settings. In addition, we present a comparison with the
existing multivariate spline methods in \cite{ALW06} and \cite{LW17} to show
that the new method produces a similar and sometimes more accurate
approximation in a more efficient fashion.
