
% Determine the 5 largest eigenvalues of K_+

format long;
type = 'pluss';    % change to type = 'minus' for K_-
lambda = 0;
opts.issym = 0;
    
L = 15;
N = 60;
betas = [2/3, 0.7:0.05:1];

lambdas = zeros(5,length(betas));
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;
    
k = 0;
for beta = betas
    k = k+1
    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);
    lambdas(:,k) = flipud(sort(diag(D)));
end

figure(5); clf;
semilogy(betas,lambdas(1,:),'k-.'); hold on;
semilogy(betas,lambdas(2,:),'k--');
semilogy(betas,lambdas(5,:),'k-');
semilogy(betas,ones(size(betas)),'k:');
v = [0.666 1 0.8 8];
axis(v);
% obvious modifications if type == 'minus'
