function [x, y, u] = sorpoisson(n,f,g,omega,tol,numiter)
%SORPOISSON Numerically approximates the solution of the Poisson
%equation on the square 0 <= x,y <= 1
%
%	[x y u] = sorpoisson(n,f,g,omega,numiter) computes the solution.
%	n is the number of subintervals, f is the right-hand side, g is
%	the boundary condition on the square, omega is the relaxation parameter,
%	and numiter is the number of SOR iterations to execute. x and y are
%	the grid of points on the x and y axes, and u is the matrix containing
%	the numerical solution. Upon building x, y, and u, sorpoisson draws a
%	surface plot of the solution.
% 
h = 1/n;
u = zeros(n+1,n+1);
x = (0:h:1)';
y = (0:h:1)';

% Assign the boundary values
u(1,:) = feval(g,zeros(n+1,1),y);
u(n+1,:) = feval(g,ones(n+1,1),y);
u(:,1) = feval(g,x,zeros(n+1,1));
u(:,n+1) = feval(g,x,ones(n+1,1));

fxy = zeros(n,n);
% Assign values of f(x,y) in the interior
for i = 2:n
	for j=2:n
		fxy(i,j) = feval(f,x(i), y(j));
	end
end

uprev = u;
k = 1;
relerr = inf;
% We have assigned the boundary values and those for f
% Now begin the SOR iteration.
while (relerr >= tol) && (k <= numiter)
	for i = 2:n
		for j = 2:n
			u(i,j) = (1-omega)*u(i,j) + ...
				omega*(u(i,j-1)+u(i,j+1)+u(i-1,j)+u(i+1,j)+h^2*fxy(i,j))/4;
		end
   end
   relerr = norm(u-uprev)/norm(u);
   uprev = u;
   k = k + 1;
end

% Plot the result
surf(x,y,u);
xlabel('x');
ylabel('y');
zlabel('z');
% use this if you want the fancy title
% title('$-\frac{\partial^{2}u}{\partial x^{2}}-\frac{\partial^{2}u}{\partial y^{2}}=5\pi^{2}\cos\left(\pi x\right)\sin\left(2\pi y\right),\\u\left(x,\, y\right)=\cos\left(\pi x\right)\,\sin\left(2\pi y\right)\,\mathrm{on\,\partial R}$',...
%       'interpreter','latex','FontSize',16,'FontWeight','bold');
