m=20;                              % number of gridpoints
A=Laplacian(m,2);                  % five point Laplacian
f=zeros(m*m,1);                    % construct right hand side
f(m:m:end)=1;
u=A\f;                             % solve by sparse Gaussian elimination
U=zeros(m+2);                      % fill solution into 2d 
U(2:m+1,2:m+1)=reshape(u,m,m);     % fills matrix column-wise
U(end,1:m+2)=1;                    % add boundary condition
x=0:1/(m+1):1; y=x;                % mesh point vectors
mesh(x,y,U);
xlabel('x'); ylabel('y');          % labeling
