function [CM_periodic] = CoeffMatrixPeriodic(K,S,PWL,alpha_nodes,mesh,QuadTriangle)
%% Construct the Stiffness matrices related to the periodic domain:

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

[Diffusion,Advection,Mass] = GlobalStiffnessMatricesPeriodic(PWL,S,mesh,QuadTriangle);

LeftBoundary  = mesh.bnd.Left;
RightBoundary = mesh.bnd.Right;
LowerBoundary = mesh.bnd.Lower;
UpperBoundary = mesh.bnd.Upper;

parfor j = 1 : alpha_number

    % Coefficient matrix
    CM_periodic{j} = - Diffusion ...
                     - 2i * Advection * alpha_nodes(j) ...
                     + Mass * (K^2 - (alpha_nodes(j))^2);

    %% Impose the periodic bnd conditions on the left and right boundaries

    CM_periodic{j}(LeftBoundary,:)...
        = CM_periodic{j}(LeftBoundary,:) + CM_periodic{j}(RightBoundary,:);

    CM_periodic{j}(RightBoundary,:) = 0;

    CM_periodic{j} = CM_periodic{j} ...
        + sparse(RightBoundary,LeftBoundary,ones(size(RightBoundary)),ng,ng) ...
        - sparse(RightBoundary,RightBoundary,ones(size(RightBoundary)),ng,ng);

    %% Impose the Dirichlet bnd conditions on bnd.Upper and bnd.Lower boundaries:

    CM_periodic{j}(LowerBoundary,:) = 0;
    CM_periodic{j}(UpperBoundary,:) = 0;

    CM_periodic{j} = CM_periodic{j} ...
        + sparse(UpperBoundary,UpperBoundary,ones(size(UpperBoundary)),ng,ng)...
        + sparse(LowerBoundary,LowerBoundary,ones(size(LowerBoundary)),ng,ng);


end
end