
% main.m - Computation of the fifth largest eigenvalue of K_+(lambda=0),
% for beta close to beta*

format long;
type = 'pluss';
lambda = 0;
opts.issym = 0;
    
Ls = 25;
betas = .91395875:2.5e-7:.91395925;

lambda5 = zeros(length(Ls),1,length(betas));
acc_checks = lambda5;
i = 0;
for L = Ls,
    i = i+1;
    Ns = L*8;
    j = 0;
    for N = Ns
        j = j+1;
        k = 0;
        h = L/N;
        if mod(N,2),
            absc = -(L/2-h/2):h:(L/2-h/2);
        else
            absc = -L/2:h:(L/2-h);
        end
        [X,Y,Z] = meshgrid(absc,absc,absc);
        R = sqrt(X.^2 + Y.^2 + Z.^2);
clear X Y Z;
        disp('Computation of the kernel of the Laplacian ...');
        % For mysinintegral.m, make sure matlab is correctly linked to
        % Maple's mfun.m
        clear Lkernel; global Lkernel; Lkernel = (1./(R+1e-20) * 1/(2*pi^2)) .* mysinintegral(pi*N*R/L);
        % The above is a higher-order approx than 1./(R+1e-20) * 1/(4*pi)
clear R;
        
        for beta = betas
            k = k+1;
            disp([num2str(i),' ',num2str(j),' ',num2str(k)]);
            global soliton;
            if beta == betas(1),
                 soliton = Petviashvili(beta, N, L);
            else
                 soliton = Petviashvili(beta, N, L, soliton);
            end
            [V,D] = eigs('BirmanSchwinger',N^3,5,'LM',opts,type,beta,lambda,L);
            d = flipud(sort(diag(D)));
            lambda5(i,j,k) = d(5);
        end
    end
end

figure(2); clf;

plot(betas-0.913959,lambda5-1)
     hold on;
plot(betas-0.913959,[0 0 0],'r-')

