clear all
close all
nmax = 8;
lambda = zeros(1,nmax);
L = 1;
% Critical value of h is -1/L
h = 5;

% Nonnegative eigenvalues
fun = @(x) tan(x)+x/(h*L);  % where x = sqrt(lambda)*L
xx = linspace(0,(2*nmax+1)*pi/2,500);
hold on
yy = tan(xx);
jump = find(abs(yy)>5);
yy(jump) = NaN;
plot(xx,yy,'r','LineWidth',2)
text(0.9*(2*nmax+1)*pi/2,-.5,'$$\sqrt{\lambda} L$$','interpreter','latex','fontsize',14)
set(gca,'XGrid','on')
set(gca,'YGrid','on')
axis([0  (2*nmax+1)*pi/2  -5  5])
plot(xx,-xx/(h*L),'r','LineWidth',2)
% Let's compute the eigenvalues numerically
disp(sprintf(' n  lambda_n^2  (2n-1)pi/2'))
for n=1:nmax
    lambda(n) = fzero(fun,[(2*n-1)*pi/2+100*eps (2*n+1)*pi/2-100*eps]);
    disp(sprintf('%2d  %10f   %10f',n, lambda(n), (2*n-1)*pi/2))
end
plot(lambda,-lambda/(h*L),'bo','LineWidth',2)
% The following is only for a "pretty" annotation of the x-axis
switch nmax
    case 1
        [hx,hy] = format_ticks(gca,...
             {'$\frac{\pi}{2}$','$\frac{3\pi}{2}$'},{'$0$'},...
             [pi/2,3*pi/2],[0],0,0,[]  ,...
             'FontSize',16,'FontWeight','Bold');
    case 2
        [hx,hy] = format_ticks(gca,...
             {'$\frac{\pi}{2}$','$\frac{3\pi}{2}$','$\frac{5\pi}{2}$'},{'$0$'},...
             [pi/2,3*pi/2,5*pi/2],[0],0,0,[]  ,...
             'FontSize',16,'FontWeight','Bold');
    case 3
        [hx,hy] = format_ticks(gca,...
             {'$\frac{\pi}{2}$','$\frac{3\pi}{2}$','$\frac{5\pi}{2}$','$\frac{7\pi}{2}$'},...
             {'$0$'},[pi/2,3*pi/2,5*pi/2,7*pi/2],[0],0,0,[]  ,...
             'FontSize',16,'FontWeight','Bold');
end
hold off

pause
% Negative eigenvalues
figure
L = 1;
h = [-.5 -1 -2 2];
xx = linspace(0,3,500);
hold on
plot(xx,tanh(xx),'r','LineWidth',2)
text(2.5,-.2,'$$\sqrt{-\lambda} L$$','interpreter','latex','fontsize',14)
axis([0  3  -2  2])
set(gca,'YGrid','on')
set(gca,'YTick',-2:2:2)
set(gca,'YTickLabel',{'-2','0','2'})
plot(xx,-xx/(h(1)*L),'r-.','LineWidth',2)
plot(xx,-xx/(h(2)*L),'r-.','LineWidth',2)
plot(xx,-xx/(h(3)*L),'r-.','LineWidth',2)
plot(xx,-xx/(h(4)*L),'r-.','LineWidth',2)
% Only h < -1/L produces a negative eigenvalue which we now compute
% numerically
h = -2;
fun = @(x) tanh(x)+x/(h*L);  % where x = sqrt(-lambda)*L, so lambda = -(x/L)^2
x_0 = fzero(fun,[0.1 3]);
disp(sprintf('lambda_0 = %10f', -(x_0/L)^2))
plot(x_0,-x_0/(h*L),'bo','LineWidth',2)
hold off