function [CM_perturbed] = B_Perturbed(K,S,PML,PWL,alpha_nodes,mesh,curve,QuadTriangle)

ng = mesh.numberpoints;
alpha_number = length(alpha_nodes);

%% Compute the coefficent of At-1 & ct-1:

[coeffDiff,coeffMass] = GradientTransformation(PML.H,QuadTriangle,curve);

%% Compute the locally perturbed periodic matrix

CM_perturbed = cell(1,alpha_number);

diffusion_perturbed = cell(3,3);
advection_perturbed = cell(3,3);
mass_perturbed      = cell(3,3);

Area_Matrix = ones(QuadTriangle.number,1) * mesh.Area;

for l = 1:3
    for m = 1:3
        
        [diff,adv,mass] = LocalStiffnessMatricesPerturbed(S,PWL,l,m,coeffDiff,coeffMass,QuadTriangle.number);
        
        diffusion_perturbed{l,m} = diff .* Area_Matrix;
        advection_perturbed{l,m} = adv  .* Area_Matrix;
        mass_perturbed{l,m}      = mass .* Area_Matrix;
        
    end
end


parfor j = 1:alpha_number
    
    Mass_perturbed = sparse([],[],[],ng,ng);
    Adv_perturbed  = sparse([],[],[],ng,ng);
    Diff_perturbed = sparse([],[],[],ng,ng);
    
   exponential = exp(conj( 1i* alpha_nodes(j)) * QuadTriangle.xnodes );

    for l = 1:3
        for m = 1:3
            
            Ig = mesh.triangles(l,:);
            Jg = mesh.triangles(m,:);
 
            % === Construct the elemenet mass and diffusion matrices
            
            element_Diff = QuadTriangle.weights *  ...
                (diffusion_perturbed{l,m} .* exponential);
            
            element_Adv = QuadTriangle.weights * ...
                (advection_perturbed{l,m} .* exponential );
            
            element_Mass= QuadTriangle.weights * ...
                (mass_perturbed{l,m} .* exponential );
            
            % === Assemble the elemenet mass and diffusion matrices
            
            Diff_perturbed = Diff_perturbed + sparse(Ig,Jg,element_Diff,ng,ng);
            Adv_perturbed  = Adv_perturbed  + sparse(Ig,Jg,element_Adv,ng,ng);
            Mass_perturbed = Mass_perturbed + sparse(Ig,Jg,element_Mass,ng,ng);
            
        end
    end
    
     CM_perturbed{j} = ( Diff_perturbed ...
        + conj(1i * alpha_nodes(j)).* Adv_perturbed - K^2 * Mass_perturbed );
    
    
    
    CM_perturbed{j}(mesh.bnd.Left,:)  = 0;
    CM_perturbed{j}(mesh.bnd.Right,:) = 0;
    
    CM_perturbed{j}(mesh.bnd.Upper,:) = 0;
    CM_perturbed{j}(mesh.bnd.Lower,:) = 0;
end