Fitting a polynomial through a given point (without lsqlin)

Has anyone written a routine to fit a polynomial constrained to run through a single arbitrary point (x0,y0) without using lsqlin from the optimisation toolbox as mentioned in http://www.mathworks.com/support/solutions/en/data/1-12BBUC/?

It is easy to "hack" the vandermonde matrix to get the polynomial to go through the origin, but it is messy to get it to go through an arbitray point. Defining new variables t = x-x0 and z = y-y0 and then fitting in the transformed variables is OK, but to recover the coefficients in x requires the pascal matrix and is very messy (but still linear). Before I do this, has anyone already attempted it?