The 2D problem of shape-preserving splines is formulated as the Differential Multipoint Boundary Value Problem (DMBVP) for thin plate tension splines. For a numerical treatment of this problem, we replace the differential operator by its difference approximation. This gives us a system of linear equations with the matrix of a special structure. We found that this matrix is positive definite. Therefore, we can solve efficiently this system of linear equations by direct or iterative methods. For the required memory of this algorithm is about n and the computational time especially the operation count is about  for each iteration, where n is the number of unknowns in the interpolation problem.