Parametric surface, normal vectors questions

Hi to all,

I have a 3D parametric defined surface that has the equation:

x=x(u,v)=(u+v)i+(u-v)j+(u^2+v^2)k

or in other way:
x=(u+v)
y=(u-v)
z=(u^2+v^2)

I have been able to draw it using only GL_POINTS and two for loops for u and v in [0,1].
Something like this:
glBegin(GL_POINTS);
for (double v=0; v <= 1; v=v+0.01 )
for (double u=0; u <= 1; u=u+0.01 )
{
{
glVertex3d(u+v,u-v,(uu+vv));
}
}
glEnd();
No problem so far (except the computational time which is very high deliberately) …
Anyway, I don’t want to use any of the primitives that OPENGL provides.
My problem is that I want to display on that surface the normal vectors.

The theory of this is:
find the cross product of the tangents (a vector) on both u,v,(first derivatives xu,xv) divide it with the length and you got the normalized normal vector on every single point of that surface.
Can someone help me with some code or links that have examples?

I don’t want to use control point defined surfaces like bezier or NURBS which I am familiar with. I just want to test ,in general, parametric defined surfaces in order to test my own algorithms for normals calculation and transformations.

Thank you in advance for any kind of help,

/Sat