function A = sysmatrix_poisson()
% function A = sysmatrix_poisson(GRID)
% fields of GRID:
% nx, ny: dimensionality
% hx, hy: mesh width
% bdry: indices of the boundary points (scalar index)
% output: system matrix A
% with ones corresponding to Dirichlet b.c. (scaled with 1/hx/hy)
% 14.08.03 by Vita
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global GRID
N = GRID.nx*GRID.ny;

b = ones(N,1);
b(GRID.bdry) = 0; % assign zeros for rows with b.c.

bu = b([N,1:N-1]);
bl = b([2:N,1]);
Dxx = spdiags([bu,-2*b,bl],[1    , 0, -1    ],N,N)/GRID.hx/GRID.hx;
Dxx = sparse(Dxx);

bu = b([N-GRID.nx+1:N,1:N-GRID.nx]);
bl = b([GRID.nx+1:N,1:GRID.nx]);
Dyy = spdiags([bu,-2*b,bl],[GRID.nx   , 0, -GRID.nx   ],N,N)/GRID.hy/GRID.hy;
Dyy = sparse(Dyy);



% boundary points
b = zeros(N,1);
b(GRID.bdry) = -1/GRID.hx/GRID.hy; % to keep the sign on diagonal
B = spdiags(b,0,N,N);

A = Dxx+Dyy+B;