S=stencil_mat(M,typ)

function S=stencil_mat(M,typ)

%this function returns the stencil matrix associated with the related
%galerkin basis set

%   1) Dirichlet/Dirichlet Mason Galerkin
%   2) Dirichlet/Dirichlet Julien Galerkin
%   3) Dirichlet/Dirichlet Trefethen Galerkin
%   4) u=u'=0
%   5) u=u''=0
%   6) Neumann/Neumann
%   7) u'[-1]-u[-1]/-1=u'[1]-u[1]==0  %stress free angular velocity
%   8) u[1]=u[-1]=u''[1]+u'[1]-u[1]=u''[-1]-u'[-1]-u[-1]=0 % stress free velocity potential
%   9) u'[1]=u'[-1]=u'''[1]=u'''[-1]==0 %Neumann/Neumann and third deriv

%create shift matrices
      L_0=sparse(eye(M));
      L_m2=sparse(diag(ones(M-2,1),-2));
      L_m4=sparse(diag(ones(M-4,1),-4));
      L_p2=sparse(diag(ones(M-2,1),2));
      
switch typ
    case 1 %Mason
        S=diag(ones(M,1));
        for j = 3:2:M-1
           S(1,j)=-1;
           S(2,j+1)=-1;
        end
        if (mod(length(M),2)==0)
            S(2,end)=-1;
        else
            S(1,end)=-1;
        end
    case 2  %Julien
        S=L_m2-L_0;
    case 3  %Trefethen
         S=diag(ones(length(g)+2,1))-diag(2*ones(length(g),1),2)+diag(ones(length(g)-2,1),4);
         S(1,3)=-1;  S(2,4)=-1;% Correct for TO and T1 mode
    case 4 %u[1]=u[-1]=u'[1]=u'[-1]=0
        D=ones(M,1);
        L2=[];
           for k=0:M-1;
              L2=[L2,-2*(k+2)/(k+3)]; 
           end
        L4=[];
           for k=0:M-1;
              L4=[L4,(k+1)/(k+3)]; 
           end 
        S=L_0*sparse(diag(D))+L_m2*sparse(diag(L2))+L_m4*sparse(diag(L4));
    case 5 %u[1]=u[-1]=u''[1]=u''[-1]=0
        D=ones(M,1);
        L2=[];
           for k=0:M-1;
              L2=[L2,-(1+((k+1)*(2*k^2+4*k+3)/((k+3)*(2*k^2+12*k+19))))]; 
           end
        L4=[];
           for k=0:M-1;
              L4=[L4,((k+1)*(2*k^2+4*k+3)/((k+3)*(2*k^2+12*k+19)))]; 
           end 
        S=L_0*sparse(diag(D))+L_m2*sparse(diag(L2))+L_m4*sparse(diag(L4));
     case 6 %u'[1]=u'[-1]==0
        D=[1];
           for k=1:M-1;
              D=[D,(k+2)^2/(2*k^2+4*k+4)]; 
           end
        L2=[0];
           for k=1:M-1;
              L2=[L2,-((k)^2)/(2*k^2+4*k+4)]; 
           end 
        S=L_0*sparse(diag(D))+L_m2*sparse(diag(L2));
     case 7 %u'[-1]-u[-1]/-1=u'[1]-u[1]==0
        D=[];
           for k=0:M-1;
              D=[D,((k+2)^2-1)/(2*(k+1)^2)]; 
           end
        L2=[];
           for k=0:M-1;
              L2=[L2,-((k^2-1))/(2*(k+1)^2)]; 
           end            
        S=L_0*sparse(diag(D))+L_m2*sparse(diag(L2));
     case 8 %u[1]=u[-1]=u''[1]+u'[1]-u[1]=u''[-1]-u'[-1]-u[-1]=0
        D=ones(M,1);
        L2=[];
           for k=0:M-1;
              L2=[L2,-2*(k+2)*(k^2+4*k+9)/((k+3)*(k^2+6*k+11))]; 
           end
        L4=[];
           for k=0:M-1;
              L4=[L4,(k+1)*(k^2+2*k+3)/((k+3)*(k^2+6*k+11))]; 
           end 
        S=L_0*sparse(diag(D))+L_m2*sparse(diag(L2))+L_m4*sparse(diag(L4));
     case 9 %u'[1]=u'[-1]=u'''[1]=u'''[-1]==0 %Neumann/Neumann and third deriv
           %%%% This stencil fails for the eigen value problem
        D=ones(M,1);
        L2=[];
           for j=0:M-1;
              L2=[L2,(-2*j^2*(11 + 2*j*(4 + j)))/((2 + j)*(3 + j)*...
                 (15 + 2*j*(6 + j)))]; 
           end 
        L4=[];
           for j=0:M-1;
              L4=[L4,(j^2*(1 + j)*(-1 + 2*j*(2 + j)))/ ...
                  ((3 + j)*(4 + j)^2*(15 + 2*j*(6 + j)))]; 
           end 
        S=L_0*sparse(diag(D))+L_m2*sparse(diag(L2))+L_m4*sparse(diag(L4));
end