figure(1)
N = 128;
[x,D] = chebdif(N,1);
I = eye(N);
while(1)
p = input('@(x) p (>0): ');
q = input('@(x) q: ');
s = input('@(x) sigma (>0): ');
A = D*(diag(p(x))*D) + diag(q(x));
b = input('@(x) left BC b(1) u + b(2) u'' = 0: ');
A(end,:) = b(1)*I(end,:) + b(2)*D(end,:);
b = input('@(x) right BC b(3) u + b(4) u'' = 0: ');
A(1,:) = b(1)*I(1,:) + b(2)*D(1,:);
B = -diag(s(x));
B(1,:) = 0;
B(N,:) = 0; 
[V,E] = eig(A,B);
E = diag(E); J = isfinite(E); E = E(J); V = V(:,J);
[F,J] = sort(abs(E)); E = E(J); V = V(:,J);
E(1:8)
for n = 1:8
    plot(x,V(:,n))
    title(['\lambda = ',num2str(E(n)/pi^2),' \pi^2'])
    E(n)/(pi*n)^2
    pause
end
end
