for k=20:1:100
    % anylitical solution
    x=(0:4/k:4); y=(0:4/k:4);
for i=1:length(x); j=1:length(y);
    u_a(i,j)=exp(y(j))*cos(x(i))-exp(x(i))*cos(y(j));
end
x0 = 0; xf = 4; % minimal and maximal value of x-coordinate
Mx = k; % number of intervals in x direction
y0 = 0; yf = 4; % minimal and maximal value of y-coordinate
My = k; % number of intervals in y direction

bx0 = inline('exp(y) - cos(y)','y'); % Boundary condition for left edge
bxf = inline('exp(y)*cos(4) - exp(4)*cos(y)','y'); % Boundary condition for right edge
by0 = inline('cos(x) - exp(x)','x'); % Boundary condition for bottom edge
byf = inline('exp(4)*cos(x) - exp(x)*cos(4)','x'); % Boundary condition for top edge

D = [x0 xf y0 yf];

MaxIter = 50000; % Maximal number of iteration
tol = 1e-6; % tollerance of Jacobi method

% Solving Laplace equation via function sollap
[u_n,x,y] = sollap(bx0,bxf,by0,byf,D,Mx,My,tol,MaxIter);

%RMS Error
U=-u_n-u_a;
U_m=mean(U);
E(k)=rms(U_m);
h(k)=4/k;
end

plot(h,E,'*')
xlabel('Grid Spacing (h)')
ylabel('Error')
grid on


function [u_n,x,y] = sollap(bx0,bxf,by0,byf,D,Mx,My,tol,MaxIter)
x0 = D(1); xf = D(2); y0 = D(3); yf = D(4); % end points for each coordinate
dx = (xf - x0)/Mx; x = x0 + [0:Mx]*dx; % step size and grid for x-coordinate
dy = (yf - y0)/My; y = y0 + [0:My]'*dy; % step size and grid for y-coordinate
Mx1 = Mx + 1; My1 = My + 1; % number of nodes for x and y directions
%boundary conditions
for m = 1:My1, u0(m,[1 Mx1])=[bx0(y(m)) bxf(y(m))]; end % left and right sides
for n = 1:Mx1, u0([1 My1],n) = [by0(x(n)); byf(x(n))]; end % bottom and top sides
u_n=u0;
sum_of_bv = sum(sum([u0(2:My,[1 Mx1]) u0([1 My1],2:Mx)']));
u0(2:My,2:Mx) = sum_of_bv/(2*(Mx + My - 2));
dx2 = dx*dx; dy2 = dy*dy; dxy2 = 2*(dx2 + dy2);
rx = dx2/dxy2; ry = dy2/dxy2;
    for itr = 1:MaxIter
    for j = 2:Mx
    for i = 2:My
u_n(i,j) = ry*(u0(i,j + 1)+u0(i,j - 1)) + rx*(u0(i + 1,j)+u0(i - 1,j)) ; % Iterative equation
    end
    end
if itr > 1 & max(max(abs(u_n - u0))) < tol, break; end % stopping condition

u0 = u_n;
    end
end