function [U,X,Y]=fdPoisson(f,g,M)

% function [u,x,y]=fdPoisson(f,g,M)
% Aim:
%     - Delta u = f in (0,1)*(0,1)
%             u = g on the boundary
% Input:
%     f,g : functions (of two variables each, vectorized)
%     M   : number of intervals in each direction
% Output:
%     U   : matrix with values of u in geographical ordering
%     X,Y : x & y coordinates of the points in geographical ordering

h=1/M;
[X,Y]=meshgrid(linspace(0,1,M+1)); 
X=X(:); Y=Y(:);

allnodes=(1:(M+1)^2)';
bdnodes=find(X==0|Y==0|X==1|Y==1);     % find indices of boundary nodes
freenodes=setdiff(allnodes,bdnodes);   % find all other indices

U=zeros(size(allnodes));

% Setting up the augmented matrix

I = speye(M+1);
E = sparse(2:M+1,1:M,1,M+1,M+1);
D = 2*I-E-E';
matrix = kron(D,I)+kron(I,D);          % -d_{xx}-d_{yy} (in this order)

% Right-hand side

U(bdnodes)=g(X(bdnodes),Y(bdnodes));
rhs = h^2*f(X,Y)-matrix*U;      % = h^2*f(X,Y)-matrix(:,bdnodes)*U(bdnodes);

% Solution and postprocessing

U(freenodes)=matrix(freenodes,freenodes)\rhs(freenodes);
U = reshape(U,M+1,M+1);
X = reshape(X,M+1,M+1);       % undoes the previous X=X(:);
Y = reshape(Y,M+1,M+1);

retur