function [a,b]=stephenson3_2d_nonunif_bd(xg,yg,nx,ny,dx,dy,alfa,beta,tim)
%
% Dalia Fishelov and Jean-Pierre Croisille
% computes the variant of 4th order of the Stephenson matrix
% for the problem alfa Delta -beta Delta ^2.
% It computes the elements of the matrix "a", which contains the
% coefficients of the operator alfa Delta -beta Delta ^2.
% It computes the right-hand side vector "b", which contains the data
% from the boundary (the right-hand side of the equation is inserted in
% another subroutine.
%
% The modification is only on the mixed term.
global xe  ye  xw  yw  xn yn  xs  ys  fle  flw  fln  fls;
global xne yne xnw ynw xse yse xsw ysw flne flnw flse flsw;
global phig;

%alfa
%beta

dim=3*nx*ny;
m1=nx-2;
n1=ny-2;
a=zeros(dim,dim);
b=zeros(dim,1);
%a=sparse(dim,dim);
%
dx2=dx*dx;
dx2inv = 1.0d0/dx2;
%
diagt = -8.0d0;
diagt1 = 2.0d0;
diagt3 =-0.5d0*dx;
%
diagb = 72.0d0;
diagb1 = -24.0d0;
diagb2 = 6.0d0;
diagb3 = 8.0d0*dx;
diagb4 = -dx;
%
diag = (alfa*diagt-beta*diagb*dx2inv)*dx2inv;
%
diag1 = (alfa*diagt1-beta*diagb1*dx2inv)*dx2inv;
diag2 = (-beta*diagb2*dx2inv)*dx2inv;
diag3 = (alfa*diagt3-beta*diagb3*dx2inv)*dx2inv;
diag4 = (-beta*diagb4*dx2inv)*dx2inv;
%
for i=2:nx-1,
  %
  for j=2:ny-1,
    %
    if(phig(i,j)==1)
      index=i+(j-1)*(m1+2);
%
      index1=3*(index-1)+1;
      index2=index1+1;
      index3=index2+1;
%      
      indexip=i+1+(j-1)*(m1+2);
      indexim=i-1+(j-1)*(m1+2);
      indexjm=i+(j-1-1)*(m1+2);
      indexjp=i+(j+1-1)*(m1+2);
%
      index1ip=3*(indexip-1)+1;
      index1im=3*(indexim-1)+1;
      index1jm=3*(indexjm-1)+1;
      index1jp=3*(indexjp-1)+1;
%
      indeximjm=i-1+(j-1-1)*(m1+2);
      indeximjp=i-1+(j+1-1)*(m1+2);
      indexipjm=i+1+(j-1-1)*(m1+2);
      indexipjp=i+1+(j+1-1)*(m1+2);
%
      index1imjm=3*(indeximjm-1)+1;
      index1imjp=3*(indeximjp-1)+1;
      index1ipjm=3*(indexipjm-1)+1;
      index1ipjp=3*(indexipjp-1)+1;
%
      a(index1,index1)=diag;
%
      a(index1,index1ip)=diag1;
%
      a(index1,index1im)=diag1;
%
      a(index1,index1jp)=diag1;
%
      a(index1,index1jm)=diag1;
%
      a(index1,index1ipjp)=diag2;
%
      a(index1,index1imjp)=diag2;
%
      a(index1,index1imjm)=diag2;
%
      a(index1,index1ipjm)=diag2;
%
      a(index1,index1ip+1)=diag3;
%
      a(index1,index1im+1)=-diag3;
%
      a(index1,index1jp+2)= diag3;
%
      a(index1,index1jm+2)=-diag3;   
%
      a(index1,index1ipjp+1)=diag4;
%
      a(index1,index1imjp+1) = -diag4;
%
      a(index1,index1ipjm+1) = diag4;
%
      a(index1,index1imjm+1) =-diag4;
%
      a(index1,index1ipjp+2) = diag4;
%
      a(index1,index1ipjm+2) = -diag4;
%
      a(index1,index1imjp+2) = diag4;
%
      a(index1,index1imjm+2) =-diag4;
%
%
%
%      
%  INSERT THE EQUATION FOR h u_x
%  -----------------------------
%
      a(index2,index2)=1.0d0;
%
      a(index2,index1ip)=-0.75d0/dx;
%
      a(index2,index1im)=0.75d0/dx;
%
      a(index2,index1ip+1)=1.0d0/4.0d0;
%
      a(index2,index1im+1)=1.0d0/4.0d0;
%      
%
%    INSERT THE EQUATION FOR h u_y
%   -----------------------------
%
      a(index3,index3)=1.0d0;
%
      a(index3,index1jp)=-0.75d0/dx;
%
      a(index3,index1jm)=0.75d0/dx;
%
      a(index3,index1jp+2)=1.0d0/4.0d0;
%
      a(index3,index1jm+2)=1.0d0/4.0d0 ;
      
  
    end;
  end;
end;
%NEW BD
%%a=zeros(dim,dim);
for j=2:ny-1,
  for i=2:nx-1,
      
    if(phig(i,j)==3)
      index=i+(j-1)*(m1+2);
      index1=3*(index-1)+1;
      index2=index1+1;
      index3=index2+1;

      indexip=i+1+(j-1)*(m1+2);
      indexim=i-1+(j-1)*(m1+2);
      indexjm=i+(j-1-1)*(m1+2);
      indexjp=i+(j+1-1)*(m1+2);
       
%
      index1ip=3*(indexip-1)+1;
      index1im=3*(indexim-1)+1;
      index1jm=3*(indexjm-1)+1;
      index1jp=3*(indexjp-1)+1;
      
      indeximjm=i-1+(j-1-1)*(m1+2);
      indeximjp=i-1+(j+1-1)*(m1+2);
      indexipjm=i+1+(j-1-1)*(m1+2);
      indexipjp=i+1+(j+1-1)*(m1+2);
%
      index1imjm=3*(indeximjm-1)+1;
      index1imjp=3*(indeximjp-1)+1;
      index1ipjm=3*(indexipjm-1)+1;
      index1ipjp=3*(indexipjp-1)+1;
%
      beta_n=2.*beta/3.;
      
%preparing alfa psi_xx- beta psi_xxxx


    h1=xg(i)-xg(i-1);
    k1=xg(i+1)-xg(i);
    if (fle(i,j)==1 )
      xxe=xe(i,j);
      yye=ye(i,j);
      k1=xxe-xg(i);
     [psies,psies_x,psies_y,psies_lap,psib_bih,inh]=sol_exa_tim_inh(xxe,yye,tim);
    end;
    if(flw(i,j)==1 )
       xxw=xw(i,j);
       yyw=yw(i,j);
       h1=xg(i)-xxw;
      [psiws,psiws_x,psiws_y,psiws_lap,psib_bih,inh]=sol_exa_tim_inh(xxw,yyw,tim);
    end;


      hpk3_1=(h1+k1)^3;
      diagt_x=2*(h1^5+k1^5-8*(h1^2*k1^3+h1^3*k1^2)-h1*k1^4-h1^4*k1)/(h1^2*k1^2);
      diagt1p_x=2*h1*(-h1^2+8*k1^2+h1*k1)/(k1^2);
      diagt1m_x=2*k1*(-k1^2+8*h1^2+h1*k1)/(h1^2);
      diagt3p_x=2*h1*(-2*k1^2-h1*k1+h1^2)/k1;
      diagt3m_x=2*k1*(2*h1^2+h1*k1-k1^2)/h1;

  
      diagb_x = 24.0d0*(h1+k1)^3/((h1*k1)^2);
      diagb1p_x = -24.0d0*(h1+3*k1)/(k1^2);
      diagb1m_x = -24.0d0*(k1+3*h1)/(h1^2);
      diagb3p_x = 24.0d0*(h1+k1)/k1;
      diagb3m_x =-24.0d0*(h1+k1)/h1;

      diag_x=(alfa*diagt_x-beta_n*diagb_x)/hpk3_1;
      diag1p_x=(alfa*diagt1p_x-beta_n*diagb1p_x)/hpk3_1;
      diag1m_x=(alfa*diagt1m_x-beta_n*diagb1m_x)/hpk3_1;
      diag3p_x = (alfa*diagt3p_x-beta_n*diagb3p_x)/hpk3_1;
      diag3m_x = (alfa*diagt3m_x-beta_n*diagb3m_x)/hpk3_1;

%preparing alfa psi_yy- beta psi_yyyy
      
      h2=yg(j)-yg(j-1);
      k2=yg(j+1)-yg(j);
      if(fln(i,j)==1)
        xxn=xn(i,j);
        yyn=yn(i,j);
        k2=yyn-yg(j);
       [psino,psino_x,psino_y,psino_lap,psib_bih,inh]=sol_exa_tim_inh(xxn,yyn,tim);
      end;
      if(fls(i,j)==1 )
        xxs=xs(i,j);
        yys=ys(i,j);
        h2= yg(j)-yys;
       [psiso,psiso_x,psiso_y,psiws_lap,psib_bih,inh]=sol_exa_tim_inh(xxs,yys,tim);
      end;

      hpk3_2=(h2+k2)^3;
      
      diagt_y=2*(h2^5+k2^5-8*(h2^2*k2^3+h2^3*k2^2)-h2*k2^4-h2^4*k2)/(h2^2*k2^2);
      diagt1p_y=2*h2*(-h2^2+8*k2^2+h2*k2)/(k2^2);
      diagt1m_y=2*k2*(-k2^2+8*h2^2+h2*k2)/(h2^2);
      diagt3p_y=2*h2*(-2*k2^2-h2*k2+h2^2)/k2;
      diagt3m_y=2*k2*(2*h2^2+h2*k2-k2^2)/h2;
      
     
      diagb_y = 24.0d0*(h2+k2)^3/((h2*k2)^2);
      diagb1p_y = -24.0d0*(h2+3*k2)/(k2^2);
      diagb1m_y = -24.0d0*(k2+3*h2)/(h2^2);
      diagb3p_y = 24.0d0*(h2+k2)/k2;
      diagb3m_y =-24.0d0*(h2+k2)/h2;
      
    
      diag_y=(alfa*diagt_y-beta_n*diagb_y)/hpk3_2;
      diag1p_y=(alfa*diagt1p_y-beta_n*diagb1p_y)/hpk3_2;
      diag1m_y=(alfa*diagt1m_y-beta_n*diagb1m_y)/hpk3_2;
      diag3p_y = (alfa*diagt3p_y-beta_n*diagb3p_y)/hpk3_2;
      diag3m_y = (alfa*diagt3m_y-beta_n*diagb3m_y)/hpk3_2;

      
%preparing  -beta\psi_{eta eta eta eta} and  -beta psi_{xi xi xi xi}
    sqrt2=sqrt(2);
    xi_ij=(yg(j)-xg(i))/sqrt2;
    eta_ij=(xg(i)+yg(j))/sqrt2;
    
%preparing  - beta psi_eta eta eta eta

    eta_ijm=(xg(i-1)+yg(j-1))/sqrt2;
    eta_ijp=(xg(i+1)+yg(j+1))/sqrt2;
    h3=eta_ij-eta_ijm;
    k3=eta_ijp-eta_ij;
    
    if (flne(i,j)==1)
       xxne=xne(i,j);
       yyne=yne(i,j);
       eta_ijp=(xxne+yyne)/sqrt2;
       k3=eta_ijp-eta_ij;
      [psine,psine_x,psine_y,psine_lap,psib_bih,inh]=sol_exa_tim_inh(xxne,yyne,tim);
       psine_eta=(psine_y+psine_x)/sqrt2;
       
    end;
    
    if(flsw(i,j)==1)
       xxsw=xsw(i,j);
       yysw=ysw(i,j);
       eta_ijm=(xxsw+yysw)/sqrt2;
       h3=eta_ij-eta_ijm;
       [psisw,psisw_x,psisw_y,psisw_lap,psib_bih,inh]=sol_exa_tim_inh(xxsw,yysw,tim);
       psisw_eta=(psisw_y+psisw_x)/sqrt2;
       
    end;
    
    hpk3=(h3+k3)^3;
    diagb_eta =  24.0d0*(h3+k3)^3/((h3*k3)^2);
    diagb1p_eta = -24.0d0*(h3+3*k3)/(k3^2);
    diagb1m_eta = -24.0d0*(k3+3*h3)/(h3^2);
    diagb3p_eta = 24.0d0*(h3+k3)/k3;
    diagb3m_eta =-24.0d0*(h3+k3)/h3;
    
      diag_eta=(-beta_n*diagb_eta)/hpk3;
      diag1p_eta=(-beta_n*diagb1p_eta)/hpk3;
      diag1m_eta=(-beta_n*diagb1m_eta)/hpk3;
      diag3p_eta = (-beta_n*diagb3p_eta)/hpk3;
      diag3m_eta = (-beta_n*diagb3m_eta)/hpk3;

   
%preparing  - beta psi_xi xi xi xi 

    xi_ijm=(yg(j-1)-xg(i+1))/sqrt2;
    xi_ijp=(yg(j+1)-xg(i-1))/sqrt2;
    h4=xi_ij-xi_ijm;
    k4=xi_ijp-xi_ij;
    
    if (flnw(i,j)==1)
      xxnw=xnw(i,j);
      yynw=ynw(i,j);
      xi_ijp=(yynw-xxnw)/sqrt2;
      k4=xi_ijp-xi_ij;
     [psinw,psinw_x,psinw_y,psinw_lap,psib_bih,inh]=sol_exa_tim_inh(xxnw,yynw,tim);
     psinw_xi=(psinw_y-psinw_x)/sqrt2;
     
  end;
  
   if(flse(i,j)==1)
      xxse=xse(i,j);
      yyse=yse(i,j);
      xi_ijm=(yyse-xxse)/sqrt2;
      h4=xi_ij-xi_ijm;
     [psise,psise_x,psise_y,psise_lap,psib_bih,inh]=sol_exa_tim_inh(xxse,yyse,tim);
     psise_xi=(psise_y-psise_x)/sqrt2;
     
  end;
  
    hpk3=(h4+k4)^3;

    diagb_xi = 24.0d0*(h4+k4)^3/((h4*k4)^2);
    diagb1p_xi = -24.0d0*(h4+3*k4)/(k4^2);
    diagb1m_xi = -24.0d0*(k4+3*h4)/(h4^2);
    diagb3p_xi = 24.0d0*(h4+k4)/k4;
    diagb3m_xi =-24.0d0*(h4+k4)/h4;
    
    diag_xi= -beta_n*diagb_xi/hpk3;
    diag1p_xi= -beta_n*diagb1p_xi/hpk3;
    diag1m_xi= -beta_n*diagb1m_xi/hpk3;
    diag3p_xi = -beta_n*diagb3p_xi/hpk3;
    diag3m_xi = -beta_n*diagb3m_xi/hpk3;
 
%inserting in a and b

      a(index1,index1)=diag_x+diag_y+diag_eta+diag_xi;

      if((fle(i,j)==1) ) 
         b(index1)=b(index1)-diag1p_x*psies;
      else
         a(index1,index1ip)=diag1p_x;
      end;

      if((flw(i,j)==1) )
         b(index1)=b(index1)-diag1m_x*psiws;
      else
         a(index1,index1im)=diag1m_x;
      end;
      
      if(fln(i,j)==1)
         b(index1)=b(index1)-diag1p_y*psino;
      else
         a(index1,index1jp)=diag1p_y;
      end;

      if(fls(i,j)==1)
         b(index1)=b(index1)-diag1m_y*psiso;
      else
         a(index1,index1jm)=diag1m_y;
      end;

     if(flne(i,j)==1) 
         b(index1)=b(index1)-diag1p_eta*psine;
      else
         a(index1,index1ipjp)=diag1p_eta;
      end;

      if(flsw(i,j)==1)
         b(index1)=b(index1)-diag1m_eta*psisw;
      else
         a(index1,index1imjm)=diag1m_eta;
      end;

      if(flnw(i,j)==1)
         b(index1)=b(index1)-diag1p_xi*psinw;
      else
         a(index1,index1imjp)=diag1p_xi;
      end;

      if(flse(i,j)==1)
         b(index1)=b(index1)-diag1m_xi*psise;
      else
         a(index1,index1ipjm)=diag1m_xi;
      end;

      if((fle(i,j)==1) )
         b(index1)=b(index1)-diag3p_x*psies_x;
      else
         a(index1,index1ip+1)=diag3p_x;
      end;

      if(flw(i,j)==1 )
         b(index1)=b(index1)-diag3m_x*psiws_x;
      else
         a(index1,index1im+1)=diag3m_x;
      end;

      if(fln(i,j)==1)
         b(index1)=b(index1)-diag3p_y*psino_y;
      else
         a(index1,index1jp+2)=diag3p_y;
      end;

      if(fls(i,j)==1)
         b(index1)=b(index1)-diag3m_y*psiso_y;
      else
         a(index1,index1jm+2)=diag3m_y;
      end;

      if(flne(i,j)==1)
         b(index1)=b(index1)-diag3p_eta*psine_eta;
      else
         a(index1,index1ipjp+1)=a(index1,index1ipjp+1)+diag3p_eta/sqrt2;
         a(index1,index1ipjp+2)=a(index1,index1ipjp+2)+diag3p_eta/sqrt2;
      end;

      if(flsw(i,j)==1)
         b(index1)=b(index1)-diag3m_eta*psisw_eta;
      else
         a(index1,index1imjm+1)=a(index1,index1imjm+1)+diag3m_eta/sqrt2;
         a(index1,index1imjm+2)=a(index1,index1imjm+2)+diag3m_eta/sqrt2;
      end;

      if(flnw(i,j)==1)
         b(index1)=b(index1)-diag3p_xi*psinw_xi;
      else
         a(index1,index1imjp+1)=a(index1,index1imjp+1)-diag3p_xi/sqrt2;
         a(index1,index1imjp+2)=a(index1,index1imjp+2)+diag3p_xi/sqrt2;
      end;
 
      if(flse(i,j)==1)
         b(index1)=b(index1)-diag3m_xi*psise_xi;
      else
         a(index1,index1ipjm+1)=a(index1,index1ipjm+1)-diag3m_xi/sqrt2;
         a(index1,index1ipjm+2)=a(index1,index1ipjm+2)+diag3m_xi/sqrt2;
      end;

%
%      
%  INSERT THE EQUATION FOR h u_x
%  -----------------------------
%
    
      diagdx=1.0d0;
      diagdxp=(h1^3*k1^2+h1^4*k1)/(h1*k1*(h1+k1)^3);
      diagdxm=(h1^2*k1^3+h1*k1^4)/(h1*k1*(h1+k1)^3);
      diagd=(2*h1^4-2*k1^4-4*h1*k1^3+4*h1^3*k1)/(h1*k1*(h1+k1)^3);
      diagdp=-(2*h1^4+4*h1^3*k1)/(h1*k1*(h1+k1)^3);
      diagdm=(2*k1^4+4*h1*k1^3)/(h1*k1*(h1+k1)^3);
      
      a(index2,index2)=1.0d0;
      
      a(index2,index1)=diagd;

      if((fle(i,j)==1) ) 
         b(index2)=b(index2)-diagdp*psies;
      else
         a(index2,index1ip)=diagdp;
      end;      

      if((flw(i,j)==1) && (phig(i-1,j)~=0)) 
         b(index2)=b(index2)-diagdm*psiws;
      else
         a(index2,index1im)=diagdm;
      end;     

       if((fle(i,j)==1) ) 
         b(index2)=b(index2)-diagdxp*psies_x;
      else
         a(index2,index1ip+1)=diagdxp;
      end;   

      if((flw(i,j)==1)) 
         b(index2)=b(index2)-diagdxm*psiws_x;
      else
         a(index2,index1im+1)=diagdxm;
      end; 

%
%    INSERT THE EQUATION FOR h u_y
%   -----------------------------
%
      
      diagdy=1.0d0;
      diagdyp=(h2^3*k2^2+h2^4*k2)/(h2*k2*(h2+k2)^3);
      diagdym=(h2^2*k2^3+h2*k2^4)/(h2*k2*(h2+k2)^3);
      diagd=(2*h2^4-2*k2^4-4*h2*k2^3+4*h2^3*k2)/(h2*k2*(h2+k2)^3);
      diagdp=-(2*h2^4+4*h2^3*k2)/(h2*k2*(h2+k2)^3);
      diagdm=(2*k2^4+4*h2*k2^3)/(h2*k2*(h2+k2)^3);
      a(index3,index3)=1.0d0;
%
      a(index3,index1)=diagd;
      if(fln(i,j)==1)
         b(index3)=b(index3)-diagdp*psino;
      else
         a(index3,index1jp)=diagdp;
      end;      
%
      if(fls(i,j)==1)
         b(index3)=b(index3)-diagdm*psiso;
      else
         a(index3,index1jm)=diagdm;
      end; 
%
      if(fln(i,j)==1)
         b(index3)=b(index3)-diagdyp*psino_y;
      else
         a(index3,index1jp+2)=diagdyp;
      end;
%
      if(fls(i,j)==1)
         b(index3)=b(index3)-diagdym*psiso_y;
      else
         a(index3,index1jm+2)=diagdym;
      end;

    end;
  end
end


