function [LAM,N,Iter,Energy] = BEC_FEM(Max_Level,mesh)

% Use the finite element method to solve the eigenvalue problem: 
% $-\Delta u + W u +\zeta u^2 u =\lambda u $ 
% with the homogeneous boundary value condition
%  INPUT: Max_Level: the number of sequence of meshes which are produced by
%  the regular refinement from  the initial mesh
%         mesh: the input mesh (the initial mesh)
%  OUTPUT: LAM: the eigenvalue approximations,
%          N: the number of elements of the sequence of meshes
%          Iter: the nonlinear iteration times in each level of mesh
%          Energy: the energy in each level mesh
% Here, we use the simple self-consistent field ("fix-point") iteration
% without damping technique. For large $\ztea$, we suggest to use the
% damping technique to guarantee the convergence.
%
% Hehu Xie: 2015-02-29

if nargin < 2
    uniformrefine = 3;
    %--------unit square: triangle-----------------------------------------
    node = [0 0; 1 0; 1 1; 0 1; 0.5 0.5];
    elem = [ 5 1 2; 5 2 3; 5 3 4; 5 4 1];
    Dirichlet = [1 2; 2 3; 3 4; 4 1];
    Neumann=[];
    %-----Create the mesh--------------------------------------------------
    mesh=Initial_FEM_matlab(node,elem,Dirichlet,Neumann);
    for i=1:uniformrefine
        mesh = refinement(mesh);
    end
end
%--------------------------------------------------------------------------
%--------------------------------------------------------------------------
coord3 = quad_triangle(13);
order = 1;
Base = Lagrange_basis_function(order);
Stiff = Init_Matrix(1,1,Base,Base,[1 1 1], [1 1 1],@Discrete_Form_Stiff,coord3);
%--------------------------------------------------------------------------
Mass = Init_Matrix(1,1,Base,Base,[1 0 0], [1 0 0],@Discrete_Form_Mass,coord3);
%------------- Assemble matrices ------------------------------------------
A1 = Matrix_Assemble(mesh,Stiff);
M1 = Matrix_Assemble(mesh,Mass);
[A1,M1] = Poission_BoundaryTreatment(mesh,Base,@u_D,A1,M1);
%-----solve a linear eigenvalue to get the initial solution----------------
%---------do the nonlinear iteration---------------------------------------
max_It = 20;
Accuracy = 1e-6;
Low = 16;
Up = 40;
Shift = 50;
%--------------------------------------------------------------------------
[xv,lm,iresult] = sptarn(A1,M1,Low+Shift,Up+Shift);
xv = xv(:,1);
[v,ind] = max(abs(xv));
xv = xv/sign(xv(ind));
xv = xv/sqrt(xv'*M1*xv);
%----------to do the nonlinear iteration-----------------------------------
Cur_Level = 1;
Stiff = Init_Matrix(1,1,Base,Base,[1 1 1], [1 1 1],@Discrete_Form_Nonlinear_Stiff,coord3);
while Cur_Level <= Max_Level    
    %----Start the nonlinear iteration-------------------------------------
    for kk = 1 : max_It
        Stiff.Aux_FEM = Init_Aux_FEM(1,xv,Base,[1 0 0]);
        xv_old = xv;
        A1 = Matrix_Assemble(mesh,Stiff);
        [A1,M1] = Poission_BoundaryTreatment(mesh,Base,@u_D,A1,M1);
        %---------Solve the eigenvalue problem-----------------------------
        [xv,lm,iresult] = sptarn(A1,M1,Low+Shift,Up+Shift);
        lm = lm(1)
        xv = xv(:,1);
        [v,ind] = max(abs(xv));
        sgn = sign(xv(ind));
        xv = xv/sgn;
        xv = xv/sqrt(xv'*M1*xv);
        res = min(norm(xv-xv_old), norm(xv+xv_old));
        fprintf('The %d-th nonlinear iteration, Res = %f\n',kk,res);
        if res < Accuracy
            fprintf('Nonlinear iteration time: %d\n',kk);
            break;
        end %end for if
        %E(Cur_Level,1) = energy_new(mesh,Base,xv,coord3,@coeff_varphi);
    end %end for kk
    Iter(Cur_Level,1) = kk;
    LAM(Cur_Level,1) = lm;
    N(Cur_Level,1) = size(mesh.elem,1);
    %------Compute the energy--------------------------------
    Energy(Cur_Level,1) = energy_new(mesh,Base,xv,coord3,@coeff_varphi);
    Cur_Level = Cur_Level + 1;
    if Cur_Level <= Max_Level
        [mesh,xv] = refinement(mesh,xv,Base);
        M1 = Matrix_Assemble(mesh,Mass);
    end
end % end for Cur_Level
LAM = LAM - Shift;
if N(end) < 2e5
    figure(1)
    showFeFunction(mesh,xv)
end
LAM
N
Iter
Energy
% print the convergence order of the eigenvalue approximations
%---use the extrapolation method to get the "accurate" eigenvalue
LAM_Exact = (4^order*LAM(end)-LAM(end-1))/(4^order-1);
LAM_Err = abs(LAM-LAM_Exact);
figure(2)
loglog(N,LAM_Err,'-o', N,0.75*LAM_Err(2)*N(2)^(order)*N.^(-order),'-.',...
       'LineWidth',3,'markersize',14)
word = ['slope=-', num2str(order)];
legend('eigenvalue errors',word);
set(gca, 'FontSize', 14)
title('Eigenvalue errors of the finite element method','fontsize',18)
xlabel('Number of elements','fontsize',18)
ylabel('Errors','fontsize',18)

%--------------------------------------------------------------------------
%  Given conditions
%--------------------------------------------------------------------------
%-----------------Sub functions-----------------------
function v = Discrete_Form_Stiff(VL,VR,X,Aux_FEM,Matrix)
v = (50 + X(:,1).^2+X(:,2).^2).*VL(:,1).*VR(:,1) + VL(:,2).*VR(:,2) + VL(:,3).*VR(:,3);
%----the discrete form for the nonlinear problem-------------
function v = Discrete_Form_Nonlinear_Stiff(VL,VR,X,Aux_FEM,Matrix)
zeta = 1;
Shift = 50;
v = (zeta*Aux_FEM(:,1).^2 +X(:,1).^2+X(:,2).^2+Shift).*VL(:,1).*VR(:,1) + VL(:,2).*VR(:,2) + VL(:,3).*VR(:,3);
%-----------------Mass Matrix---------------------
function v = Discrete_Form_Mass(VL,VR,X,Aux_FEM,Matrix)
v = VL(:,1).*VR(:,1);
%-----------Coefficient function degree 0----------------------------------
function v = coeff_varphi(x)
%--------Model problem-----------------------------------------------------
%-------Example 1----------------------------------------------------
v = x(:,1).^2+x(:,2).^2;
%--------------------------------------------------------------------------
function v=u_D(x)
v = sparse(size(x,1),1);
%--------------------------------------------------------------------------
function v=u_N(x)
v = sparse(size(x,1),1);