In 13.4.3, Multitaper methods, equation 13.4.18 describes a trigiagonal matrix whose eigenvectors are the discrete prolate spheriodal sequences which are used as window functions for spectral estimation.

If I plug-in this equation into Mathematica, I get eigenvectors that are different than what the NR code returns. I am pretty sure I am doing something wrong, but I wanted to check that there is no typo here. The implementation in the NR code seems consistent with this equation.

Mathematica code:
s = SparseArray[{
{i_, j_} /; i == j ->
1/4 (n^2 - (n - 1 - 2 (i - 1))^2 Cos[2 Pi jres/n]),
{i_, j_} /; i == j - 1 -> -1/2 i (n - i),
{i_, j_} /; i == j + 1 -> -1/2 j (n - j)}, {n, n}] // N;

Array indexing in Mathematica starts at 1.

Gruenbacher and Hummels, IEEE Sig Proc, vol 42, no 11, 1994, gives a similar formulation for this matrix and I couldn't get this one to work either:

t = SparseArray[{
{i_, j_} /; i == j -> ((n - 1)^2 - (i - 1))^2 Cos[2 Pi jres/n],
{i_, j_} /; i == j - 1 -> 1/2 i (n - i),
{i_, j_} /; i == j + 1 -> 1/2 j (n - j)}, {n, n}] // N;