poulson

  function M = Inverter(nx,ny)
  
  A=eye(nx*ny,nx*ny);
  for x=1:nx, for y=1:ny,
    i=x+(y-1)*nx;
    if( x ~= 1 ), A(i,i-1)=1; end;
    if( x ~= nx ), A(i,i+1)=1; end;
    if( y ~= 1 ), A(i,i-nx)=1; end;
    if( y ~= ny ), A(i,i+nx)=1; end;
  end; end;

  b=ones(nx*ny,1);
  
  s=gflineq(A,b);

  M=zeros(ny,nx);
  for x=1:nx, for y=1:ny, M(x,y)=s(x+(y-1)*nx); end; end;
  
  return