I don't think those equations are right... I just tried them now in the 2D case, and they didn't give the right answer at all. (I used an ellipse with a = 2 and b = 1, and (x0, y0) = (2, 1.5), and plotted the parametric curve according to the equation you posted. It goes through a point on the ellipse, but definitely not point on the ellipse nearest to (2, 1.5).)
But, assuming that some typo was made, I believe I can give a partial explanation:
Imagine taking the point (x0, y0, z0) and drawing a small sphere around it, then gradually expanding that sphere till it just touches the ellipsoid. The point where it touches the ellipsoid is the place you're looking for, and since the sphere and ellipsoid are just touching at that point, they're tangent there. This means the ellipsoid's normal at that point is in the opposite direction to the sphere's normal - and since the sphere's normal points directly away from its center, it follows the ellipsoid's normal points directly towards the sphere center, aka (x0, y0, z0).
So we know the location we want is the point on the ellipsoid where its normal points toward (x0, y0, z0). Now, if you define an ellipsoid as a scalar field f(x,y,z) where f = 0 for points on the surface, f \< 0 inside and f > 0 outside, then the normal of an ellipsoid is given by the gradient of f.
Now, here is where I'm not sure what your author is doing anymore - he seems to be trying to construct a 1D curve that (somehow) will go through the right point on the ellipsoid. At this point it makes more sense to me to parameterize the ellipsoid as a function P(u,v) mapping u,v to some point on the ellipsoid surface. Then, what we want to find are the u,v for which
(x0, y0, z0) - P(u,v) = k * gradientF(P(u,v))
for some k, i.e. the vector from the point on the surface to (x0, y0, z0) is in the same direction as the gradient at that point. This can be re-expressed as a function that normalizes both vectors and takes their dot product:
G(u,v) = normalize((x0, y0, z0) - P(u,v)) dot normalize(gradientF(P(u,v)))
Now, this function will have a maximum exactly at the u,v coordinates where you'll find the point you're looking for. You have to apply a two-dimensional maximization now, but that's not so bad as the function is very well behaved, there are no local maxima. You could even probably solve it analytically, though it would be messy (and possibly the numerical solution would perform better).