p21.m

% p21.m - eigenvalues of Mathieu operator -u_xx + 2qcos(2x)u
% (compare p8.m and p. 724 of Abramowitz & Stegun)
N = 42; h = 2*pi/N; x = h*(1:N);
D2 = toeplitz([-pi^2/(3*h^2)-1/6 ...
-.5*(-1).^(1:N-1)./sin(h*(1:N-1)/2).^2]);
qq = 0:.2:15; data = [];
for q = qq;
e = sort(eig(-D2 + 2*q*diag(cos(2*x))))';
data = [data; e(1:11)];
end
clf, subplot(1,2,1)
set(gca,'colororder',[0 0 1],'linestyleorder','-|--'), hold on
plot(qq,data,'linewidth',.8), xlabel q, ylabel \lambda
axis([0 15 -24 32]), set(gca,'ytick',-24:4:32)