Software package for plotting 3-d splines Given a finite point set $P \subset \mathbb{R}^2$ and a height function $h:P -> \mathbb{R}$ I want to produce a smooth surface that interpolates between the values $\left\{[p~h(p)]^\top \in \mathbb{R}^3 |~p \in P\right\}$. 
Question (unrefined): What software packages are out there for computing and visualizing such surfaces from a list of points?
I know that the Sage method list_plot3d can be used with the interpolation_type=spline but this command does not work as expected. For example, inputting 

list_plot3d([[1, 1, 1], [1, 2, 1], [0, 1, 3], [1, 0, 4]], point_list=True, interpolation_type='spline')

into Sage returns the error

TypeError: m >= (kx+1)(ky+1) must hold

(A discussion about properly formatting the input into sage can be found here.)
Question (refined): What open source software exists for computing and visualizing a polynomial surface from a finite list of points in $\mathbb{R}^3$?
 A: Your points $P$ do not appear to lie on a rectangular grid. So, they would be referred to as "scattered", and what you're doing is "scattered data interpolation". If you search for this term, you'll get lots of hits. One example is this Wikipedia page.
If you want an interpolating polynomial, you simply write down a polynomial in which the number of coefficients is the same as the number of your data points. Then, for each point of $P$, you write down an equation that expresses the interpolation. You get a system of linear equations, which you can (usually) solve to get the coefficients of the polynomial. I say "usually" because sometimes the linear systems arising in 2D interpolation don't have unique solutions.
If the number of points is large, you might run into problems -- the so-called Runge phenomenon will bite you. 
If you want a software package that does the interpolation for you, then there are various recommendations in the answers listed below. Matlab seems to be people's favorite solution:
Link1
Link2
Link3
Link4
I have never used SAGE, but I can guess what the error message is telling you. A tensor product spline (or polynomial) surface of degree $kx \times ky$ will have at least $(kx+1) \times (ky+1)$ coefficients. In order for some internally constructed system of equations to be solvable (as I outlined above), the number of input data points $m$ must be at least this large. The "spline" option probably tells SAGE to use cubic splines, which means you will need at least 16 data points.
