I am aware of how to compute the cubic spline of a one-variable function with mm control points using a linear equation . However, I am not sure how to do any type of spline on a two-variable function with m knots. Is this even possible?
The following sources may lead you (though, not directly) to the algorithmsyou need
a book:
Bases in Function Spaces, Sampling, Discrepancy, Numerical Integration, by Hans Triebel
a review
https://pdfs.semanticscholar.org/e9d3/1e439cfe3430d3388b09ddf11c526c3453c2.pdf, by C. de Boor
a paper
Tensor products of Sobolev–Besov spaces and applications to approximation from the hyperbolic cross, December 2009, Journal of Approximation Theory 161(2):748-786, DOI: 10.1016/j.jat.2009.01.001 by Winfried Sickel &Tino Ullrich
a paper by G. G. Lorentz & Rudolph A. Lorentz
Article Bivariate Hermite Interpolation and Applications to Algebraic Geometry