function [A,p,b] = pois2d(x,y,M,N,condition)

size = M*N;

dx = x/M;
dy = y/N;

dx_squared = dx^2;
dy_squared = dy^2;

if (strcmp(condition, 'pure_Neumann'))
  A = spalloc(size - 1, size - 1, 5 * (size - 1));
elseif (strcmp(condition, 'partly_Dirichlet'))
  A = spalloc(size - M, size - M, 5 * (size - M));
end

p = zeros(length(A), 1);
b = dx_squared * dy_squared * rand(length(A), 1);

for i=1:size-1,

  % 1st corner node
  if (i == 1)
      A(i,i) = dx_squared + dy_squared;
      A(i, i+1) = - dy_squared;
      A(i, i+M) = - dx_squared;
  % 2nd corner
  elseif (i == M)
      A(i,i) = dx_squared + dy_squared;
      A(i, i-1) = - dy_squared;
      A(i, i+M) = - dx_squared;
  % 3rd corner
  elseif (i == size - M + 1)
    if (strcmp(condition, 'pure_Neumann'))
      A(i,i) = dx_squared + dy_squared;
      A(i, i+1) = - dy_squared;
      A(i, i-M) = - dx_squared;
    end
  % 4th corner removed; fix adjoining nodes
  elseif (i == size - M)
      A(i,i) = 2*dx_squared + dy_squared;
      A(i, i-1) = - dy_squared;
      A(i, i-M) = - dx_squared;
  elseif (i == size - 1)
    if (strcmp(condition, 'pure_Neumann'))
      A(i,i) = dx_squared + 2*dy_squared;
      A(i, i-1) = - dy_squared;
      A(i, i-M) = - dx_squared;
    end
  % fix node right below the 3rd corner
  elseif (i == size - 2*M +1)
      A(i,i) = 2*dx_squared + dy_squared;
      A(i, i+1) = - dy_squared;
      A(i, i-M) = - dx_squared;
    if (strcmp(condition, 'pure_Neumann'))
      A(i, i+M) = - dx_squared;
    end
  % lower horizontal boundary?
  elseif ((i >= 1) & (i <= M))
      A(i,i) = dx_squared + 2*dy_squared;
      A(i, i-1) = - dy_squared;
      A(i, i+1) = - dy_squared;
      A(i, i+M) = - dx_squared;
  % upper horizontal boundary?
  elseif ((i > size - M) & (i <= size))
    if (strcmp(condition, 'pure_Neumann'))
      A(i,i) = dx_squared + 2*dy_squared;
      A(i, i-1) = - dy_squared;
      A(i, i+1) = - dy_squared;
      A(i, i-M) = - dx_squared;
    end
  % nodes right below upper horizontal boundary
  elseif ((i > size - 2*M) & (i <= size - M))
      A(i,i) = 2*(dx_squared + dy_squared);
      A(i, i-1) = - dy_squared;
      A(i, i+1) = - dy_squared;
      A(i, i-M) = - dx_squared;
    if (strcmp(condition, 'pure_Neumann'))
      A(i, i+M) = - dx_squared;  
    end

  % left vertical boundary?
  elseif (rem(i-1, M) == 0)
      A(i,i) = 2*dx_squared + dy_squared;
      A(i, i+M) = - dx_squared;
      A(i, i-M) = - dx_squared;
      A(i, i+1) = - dy_squared;
  % right vertical boundary?
  elseif (rem(i, M) == 0)
      A(i,i) = 2*dx_squared + dy_squared;
      A(i, i+M) = - dx_squared;
      A(i, i-M) = - dx_squared;
      A(i, i-1) = - dy_squared;
  % internal node?
  else 
      A(i,i) = 2*dx_squared + 2*dy_squared;
      A(i,i+1) = - dy_squared;
      A(i,i-1) = - dy_squared;
      A(i, i+M) = - dx_squared;
      A(i, i-M) = - dx_squared;
  end
end