Rendering Polynomial surfaces in OpenGL

I need to be able to draw polynomial surfaces but I have
not found a way to do it with OpenGL. The problem is as
follows:

Given the polynomial: f(x,y) = sum(a_ijx^iy^j) i,j = 1,2,…,n
example: f(x,y) = 2.0x^2 + y^2 +3xy + x + 2.0y + 1.0

I want to be able to draw the 3D surface over a triangular domain:
example: I have the 2D triangle T2=(x1,y1, x2,y2, x3,y3).
If I say something like:

glBegin(GL_TRIANGLES)
glVertex3f(x1,y1,f(x1,y1));
glVertex3f(x2,y2,f(x2,y2));
glVertex3f(x3,y3,f(x3,y3));
glEnd();

I will see the the 3D triangle T3 whose vertices are
(xi,yi,f(xi,yi)) i=1,2,3 and all the points in the
2D triangle T2 have z coordinate such that it
lies inside the plane of the the triangle T3.

What I want to change is that each point P2(xp,zp) inside the
triangle T2 is mapped to a point P3(xp,yp,zp) where zp is given
by: zp = f(xp,yp).

How can I achieve that with OpenGL?

Thanks in advance!

Either 1) use smaller triangles so it more closely approximates your function

or 2) learn a whole lot about fragment programs.

I suggest 1).

Or 3) Split your triangle into many smaller ones and use a vertex program to calculate the polynomial (a lot less magic required compared to fragment programs). You would also need to calculate normals, which can be done in the vp too.

Honestly I don’t think you can do this with fragment programs unless you do the ‘ray trace in fp’ trick, which really is pretty involved.