% This Matlab code is designed for solving
%
%   -(u_xx + u_yy) + alf u_x = f(x,y),  a <= x<=b,  c<=y <= d
%
%  with Dirichlet BC at x=a, y=c, y=d, and Neumann BC at x=b using SOR.
%
%  Functions needed: f3(x,y) which is f(x,y), ux3(x,y) the Neumann BC,
%  uexact3(x,y), an exact solution for testing purpose.

clear all; close all;

global alf

n  = input('Enter n = '); m=n;
w = input('Enter w= ');
tol = input('Enter tol = ');

a =1; b=3; c = 2; d=4; alf = 1.5; h=(b-a)/n; h1=h*h;
x = a:h:b; x=x'; y=c:h:d; y=y'; 

error = 1000; u0 = zeros(m+1,n+1);

% Set boundary condition

for i=1:m+1
  u0(i,1) = uexact3(x(i),c);
  u0(i,n+1) = uexact3(x(i),d);
end

for j=1:n+1,
  u0(1,j) = uexact3(a,y(j));
%  u0(m+1,j) = uexact3(b,y(j));    % For Dirichlet BC
end

k = 0; u=u0;
while error > tol
  for i=2:m
     for j=2:n
	u_temp = (u(i-1,j) + u(i+1,j) + u(i,j-1) + u(i,j+1) )/4 ...
                 -alf*h*(u(i+1,j)-u(i-1,j))/8 + h1*f3(x(i),y(j))/4;
	u(i,j) = (1-w)*u0(i,j)  + w*u_temp;
     end
  end

  for j=2:n
     u_temp = (2*u(m,j) + u(m+1,j-1) + u(m+1,j+1) + 2*h*ux3(b,y(j)) )/4 ...
                  -alf*h1*ux3(b,y(j))/4 + h1*f3(b,y(j))/4;
         u(m+1,j) = (1-w)*u0(m+1,j)  + w*u_temp;
   end

  error = max(max(abs(u-u0)));
  u0 = u;
  k = k+1;
end

for i=1:m+1
  for j=1:n+1
    ue(i,j) = uexact3(x(i),y(j));
  end
end

figure(1); mesh(x,y,u); title('The computed solution');
figure(2); mesh(x,y,u-ue); title('The error plot');

Abs_error = max(max(abs(ue-u))), 
Rel_error = Abs_error/max(max(abs(ue))),
k