the canonical example of a Morse function is the height function on a torus

preliminaries

Needs["Graphics`ParametricPlot3D`"]

Needs["Calculus`VectorAnalysis`"]

Get["Analysis`VectorField.m`"]

Needs["Graphics`Shapes`"]

<<"/xinwei/mma/Topology/Holes.ma"

also "/xinwei/mma/Topology/Holes.ma"

toro[a_,b_] = {(2+Cos[b]) Cos[a], Sin[b], (2+Cos[b]) Sin[a] };

ParametricPlot3D[ toro[a,b], {a,0,2Pi}, {b,0,2Pi}, Axes->Automatic]

define tangent & normal spaces

nt[a_,b_] = {Cos[b] Cos[a], Sin[b], Cos[b] Sin[a] };

ta[a_,b_] = D[toro[a,b],a];
tb[a_,b_] = D[toro[a,b],b];

{nt[a,b].ta[a,b], nt[a,b].tb[a,b]}

check that nt is indeed unit length

nt[a,b]. nt[a,b]

SetCoordinates[Cartesian];

Grad[z]

is constant on R3, project onto tangent space...

Df[a_,b_] = Grad[z] - Grad[z].nt[a,b] nt[a,b]

Df[a_,b_] = {-(Cos[a]*Cos[b]^2*Sin[a]), -(Cos[b]*Sin[a]*Sin[b]), 
  1 - Cos[b]^2*Sin[a]^2};

is the tangential gradient of f at <a,b> on the torus.
In local coordinates of the torus, f is defined

the Hessian of f
in local coordinates, f can be defined

f[a_,b_] = (toro[a,b][[3]])

so we can just compute the Hessian of f with respect to local coords

MatrixForm[Hf[a_,b_] = {{D[f[a,b],{a,2}],D[f[a,b],a,b]},
				{D[f[a,b],a,b],D[f[a,b],{b,2}]}}]

Hf[a,b]

Hf[a_,b_] = {{-((2 + Cos[b])*Sin[a]), -(Cos[a]*Sin[b])}, 
  {-(Cos[a]*Sin[b]), -(Cos[b]*Sin[a])}};

hessf[a_,b_] = Det[Hf[a,b]]

hessf[a_,b_] = 2*Cos[b]*Sin[a]^2 + Cos[b]^2*Sin[a]^2 - Cos[a]^2*Sin[b]^2;

We can cut this to RealTimeAlgebra (in NexTDeveloper/Demos)

SolveAlways[hessf[a,b] == 0, {a,b}]

Look at critical points of f, where the gradient f vanishes:

SolveAlways[Df[a,b] == {0,0,0},{a,b}]

{{Cos[b]^2 -> Sin[a]^(-2), Cos[a] -> 0, Sin[b] -> 0}, 
  {Cos[b]^2 -> Sin[a]^(-2), Cos[a] -> 0, Sin[b] -> 0}}

toro[a,b] /.%30[[1]]

There are four solutions , thus four critical points on the torus

MatrixForm[%35 //.
{{Cos[b]->1,Sin[a]->1},{Cos[b]->-1,Sin[a]->1},
{Cos[b]->1,Sin[a]->-1},{Cos[b]->-1,Sin[a]->-1}}]

A critical point is nondegenerate if the hessian is  nondegenerate at that point.

Let's look again at the torus

tor = ParametricPlot3D[ toro[a,b], {a,0,2Pi, Pi/12}, {b,0,2Pi, Pi/12}]

and take slices of this torus by level sets of  f, in this case simply the planes z= constant.
(we'll use  the ExtrinsicCut op from "Holes")

vol[z_] = {{-3,3},{-1,1},{z,4}};

Do[
	ExtrinsicCut[tor,vol[zz]],
	{zz,2.9,0,-.2}
	] //Timing

show only two frames; even at this coarse resolution, the pictures occupy too much memory.  (Currently, the NeXT version does not support graphics conversions which are available on the Macintosh version 1.2.)

{1039.483*Second, Null} (15 frames on the NeXT)

let's look at the (tangential) vector field corresponding to the gradient of f.

VectorSurfaceAnalytic[

VectorParametricSurface[ function_List, field_List, 
		{x_, x0_, x1_, dx_:((x1-x0)/(plotpoints-1))}, 
		{y_, y0_, y1_, dy_:((y1-y0)/(plotpoints-1))},
		scaling_,
		options___] := Show[Graphics3D[
		Table[vectorLinesSurface[function, 
					field 
					scaling],
				{x, x0, x1, dx}, {y, y0, y1, dy}]
	  ], options]

VectorParametricSurface[toro[x,y],Df[x,y],  {a,0,2Pi, Pi/12}, {b,0,2Pi, Pi/12},
	0.8, Boxed->False,PlotRange->{{-3,3},{-1,1},{-3,3}}]