Based on the image you provided, I would say that you are on the right track, and the next step is to enable smooth shading on the mesh. Currently, you are seeing hard edges in the lighting on adjacent triangles because each triangle is being flat shaded. The normal of only one of the vertices is being used to do lighting calculations, rather than the normals of all three vertices being interpolated across the face of the triangle (or the lighting at each vertex being interpolated, depending on whether you're using per-pixel or per-vertex lighting, respectively).
In older OpenGL implementations, enabling smooth (Gourard) shading would be done with the following line:
but note that this is now deprecated in favor of the programmable pipeline. Fear not, though, as it will still work if you work with an older OpenGL library (and it may also be supported in the most recent version; someone may correct me if I'm wrong about that). To go off on a brief tangent, learning OpenGL via the deprecated immediate mode is sometimes preferred because it's a little more straightforward than the programmable pipeline.
Anyway, the smooth shading is also dependent upon the normals of your vertices being set correctly, which seems to already be happening, upon examination of your image. If, for whatever reason, they are not quite correct, I'd be happy to update this answer to provide some direction on how to calculate them.
Hope this helps :)
Update: To calculate the normals for proper lighting, you can use the following rule: the normal for a vertex is equal to the average of the face normals of the triangles it is a part of. The face normal of a triangle is calculated by taking the cross-product of two of its edges. So, for one of the corners of your mesh:
| .` |
You can calculate the face normal of
ABC by taking
BA x BC, and the face normal of
BCD to by taking
BC x BD. Averaging the two resulting vectors together and normalizing the result gets you the normal for
B. For vertices in the middle of the mesh, you will be calculating and averaging the face normals of eight triangles. Note that this is the method which bobobobo provided in his answer. I hope that this summary acts as a good complement to it.
For this kind of calculation, it really pays to store your mesh in a Half-Edge data structure. This makes traversal of the triangles much easier, and also allows you to implement subdivision and decimation of the mesh in case you need it later on.