function [LAM,LAM_Dir,Dis_1,Dis_0,Err,Err_dir,N]=Full_Eigen_Multigrid(Max_Level,mesh)

tic;
if nargin < 2
    uniformrefine = 4;
    %--------unit square: triangle---------------------------------------------
    node = [0 0; 1 0; 1 1; 0 1];
    elem = [ 1 2 4;
             2 3 4];  
    Dirichlet = [1 2;2 3; 3 4; 4 1];
    Neumann=[];
    %--------reference triangle------------------------------------------------
    % node =[0 0;1 0;0 1];
    % elem=[1 2 3 -1];        
    % Dirichlet=[1 2; 2 3; 3 1];
    % Neumann=[];
    mesh=Initial_FEM_matlab(node,elem,Dirichlet,Neumann);
    for i=1:uniformrefine
        mesh = refinement(mesh);
    end
end

%---------First we refine the mesh to build the multilevel mesh------------
coord3 = quad_triangle(13);
order = 1;
Base = Lagrange_basis_function(order);
%-----solve the eigenvalue in the linear finite element--------------------
Stiff = Init_Matrix(1,1,Base,Base,[0 1 1], [0 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,bd] = Poission_BoundaryTreatment(mesh,Base,@u_D,A1,M1);
%-----eigenvalue-----------------------------------------------------------
[xv,lm,iresult] = sptarn(A1,M1,0,100);
%-----------to do the correction first time--------------------------------
lm = lm(1);
xv = xv(:,1);
[v_max,index] = max(abs(xv(:,1)));
v_max = xv(index,1);
xv = xv/v_max;
LAM(1) = lm;
LAM_Dir(1) = lm;

Dis_1(1) = 0;
Dis_0(1) = 0;
Err(1) = error_estimate(mesh,Base,coord3,xv,@Exact_Grad,'gradient');
Err_dir(1) = Err(1);
%-----Start the full multigrid method--------------
mg = MGInit(A1,M1,bd);
N(1) = size(mesh.elem,1);
%-----------Save the initial mesh------------------------------------------
%-----------to do the multilevel correction method-----
mesh = Init_Multi_Mesh(mesh);
meshs = [mesh,mesh];
clear mesh;
Cur_Level = 2;
%---------Multigrid eigenvalue solution------------------------------------
Num_Correction = 2*ones(Max_Level,1);
Num_Correction(end) = 2;
Num_Eig = 1;

while Cur_Level <= Max_Level
    Cur_Level
    %------------propolongation---------------------`
    [meshs(1),xv,Prolong] = refinement(meshs(1),xv,Base);
    xv = xv(:,1);
    N(Cur_Level) = size(meshs(1).elem,1);
    %------First do the eigenvalue problem solving in the coarse space-----
    %-----build the stiff matrix ------------------------------------------
    %-----build the stiff matrix ------------------------------------------
    Stiff.R_Mesh = 1;
    Mass.R_Mesh = 1;
    A = Matrix_Assemble(meshs(1),Stiff);
    M = Matrix_Assemble(meshs(1),Mass);
    [A,M,bd] = Poission_BoundaryTreatment(meshs(1),Base,@u_D,A,M);
    %------Direct solving ----
    [xv_dir,lm_dir,iresult] = sptarn(A,M,19,40);
    lm_dir = lm_dir(1);
    xv_dir = xv_dir(:,1);
    [v_max,index] = max(abs(xv_dir(:,1)));
    v_max = xv_dir(index,1);
    xv_dir = xv_dir/v_max;
    mg = MGAddLevel(mg,A,M,Prolong,bd);
    Stiff.R_Mesh = 2;
    Mass.R_Mesh = 2;
    AHh = Matrix_Assemble(meshs,Stiff);
    MHh = Matrix_Assemble(meshs,Mass);
    %----------Do the correction---------------
    for iter = 1 : Num_Correction(Cur_Level)
        RHS = M*xv*diag(lm);
        [A,RHS] = Poission_BoundaryTreatment(meshs(1),Base,@u_D,A,RHS);
        xv = MGIteration(mg,xv,RHS,Cur_Level,1,2);
        %---------to do the correction second time-------------------------
        ann = xv'*A*xv;        
        RHS = AHh*xv;
        AA = [A1,RHS; RHS',ann];        
        RHS = MHh*xv;
        MM = [M1,RHS;RHS',xv'*M*xv];
        [AA,MM] = Poission_BoundaryTreatment(meshs(2),Base,@u_D,AA,MM);
        %-- solve the eigenvalue problem ----------------------------------
        [xv2,lm,iresult] = sptarn(AA,MM,18,100);        
        lm = lm(1)
        xv2 = xv2(:,1);
        iter
        xv = Get_Solution_Multilevel(Base,meshs(2),xv2(1:end-Num_Eig,:),Base,...
                                     meshs(1),xv2(end-Num_Eig+1:end,:),xv);  
        [v_max,index] = max(abs(xv(:,1)));
        v_max = xv(index,1);
        xv = xv/v_max;
        if Cur_Level == Max_Level
            %  compute the error estimates in the finest level of space
                Dis_1(iter) = difference_estimate(meshs(1),Base,xv_dir,Base,xv,...
                                         coord3,1);
                Dis_0(iter) = difference_estimate(meshs(1),Base,xv_dir,Base,xv,...
                                         coord3,0);
                Dis_LAM(iter) = abs(lm-lm_dir);
        end
   end  
    LAM(Cur_Level) = lm;
    LAM_Dir(Cur_Level) = lm_dir;
    Err(Cur_Level) = error_estimate(meshs(1),Base,coord3,xv,@Exact_Grad,'gradient');
    Err_dir(Cur_Level) = error_estimate(meshs(1),Base,coord3,xv_dir,@Exact_Grad,'gradient');                   
    Cur_Level = Cur_Level + 1;    
end
r_1 = Dis_1(end)/Dis_1(1);
r_1 = exp(log(r_1)/(Num_Correction(end)-1));
r_0 = Dis_0(end)/Dis_0(1);
r_0 = exp(log(r_0)/(Num_Correction(end)-1));
r_LAM = Dis_LAM(end)/Dis_LAM(1);
r_LAM = exp(log(r_LAM)/(Num_Correction(end)-1));
X = 1: Num_Correction(end);
semilogy(X,Dis_1, '-*', X,Dis_0,'-*', X, Dis_LAM, '-s',...
         X, Dis_1(2)*r_1^(-X(2))*r_1.^(X),'-.', X, Dis_0(2)*r_0^(-X(2))*r_0.^(X),'-.',...
         X, Dis_LAM(2)*r_LAM^(-X(2))*r_LAM.^(X),'-.','LineWidth',3,'markersize',14)
%--------------------------------------------------------------------------
word1 = ['line: ', num2str(r_1),'^x'];
word2 = ['line: ', num2str(r_0), '^x'];
word3 = ['line: ', num2str(r_LAM), '^x'];
legend('|u_h-u_h^{dir}|_1','||u_h-u_h^{dir}||_0','|\lambda_h-\lambda_h^{dir}|', word1,word2,word3);
set(gca, 'FontSize', 18)
title('Errors of the finite element method','fontsize',18)
xlabel('Number of elements','fontsize',18)
ylabel('Errors','fontsize',18)
toc

%--------------------------------------------------------------------------
%  Given conditions
%--------------------------------------------------------------------------
%-----------------Sub functions-----------------------
function v = Discrete_Form_Stiff(VL,VR,X,Aux_FEM,Matrix)
v = VL(:,1).*VR(:,1) + VL(:,2).*VR(:,2);
%-----------------Mass Matrix---------------------
function v = Discrete_Form_Mass(VL,VR,X,Aux_FEM,Matrix)
v = VL(:,1).*VR(:,1);
function v = Exact(x)
v=sin(pi*x(:,1)).*sin(pi*x(:,2));

function v = Exact_Grad(x)
v = pi*[cos(pi*x(:,1)).*sin(pi*x(:,2)),sin(pi*x(:,1)).*cos(pi*x(:,2))];

function v=coeff(x)
epsilon=1;
dim=size(x,1);
v=zeros(dim,2,2);
v(:,1,1)=epsilon;
v(:,1,2)=0;
v(:,2,1)=0;
v(:,2,2)=epsilon;

function v=mass(x)
v =ones(size(x,1),1);

function v=RHS(x)
v =  2*pi*pi*sin(pi*x(:,1)).*sin(pi*x(:,2));
%--------------------------------------------------------------------------
function v=u_D(x)
v = zeros(size(x,1),1);

function v=u_N(x)
v = zeros(size(x,1),1);