function [] = testtoeplitztridiagonal(),
   al = 1.234567; n = 5;
   T = ToeplitzTridiagonal(al, n);
   Tinv = InverseOfToeplitzTridiagonal(al, n);
   disp('T ='); display(T);
   disp('Tinv ='); display(Tinv);
   P = T * Tinv; P(abs(P) < 1.0E-15) = 0;
   disp('T * Tinv after rounding small numbers to zero = '); display(P);
   disp('norm(T * Tinv - eye(n)) = '); display(norm(T * Tinv - eye(n)));
end

function [T] = ToeplitzTridiagonal(al, n),
   col = [2 + al; -1; zeros(n - 2, 1)];
   T = toeplitz(col);
end

function [Tinv] = InverseOfToeplitzTridiagonal(al, n),
   be = 2 * log(sqrt(al) / 2 + sqrt(1 + al / 4));
   coef = 1 / (2 * sinh(be) * sinh((n + 1) * be));
   ind = (0 : n - 1)';
   coltoeplitz = coef * cosh((n + 1 - ind) * be);
   colhankel = coef * cosh((n - 1 - ind) * be);
   S = toeplitz(coltoeplitz);
   H = hankel(colhankel, colhankel(n : -1 : 1));
   Tinv = S - H;
end

