function basic_iterative_methods()
figure(1);
hold on;
p1 = 4;
for p = p1 + 1 : 10
    n = 2^p + 1;
    n1 = n - 1;
    t = linspace(0,1,n);
    [x,y] = meshgrid(t,t);
    f = 2*(x.^2 - x + y.^2 - y);
    h(p - p1) = 1/(n - 1);
    h2 = h(p - p1)^2;
    u = zeros(n);
    u_exact = x.*y.*(1 - x).*(1 - y);

    itermax = 150;
    itermin = 100;
    ii = [itermin : itermax];
    % Jacobi
    I = [2 : n1];
    Ip = [3 : n];
    Im = [1 : n - 2];
    u = zeros(n);
    unew = u;
    tic
    for k = 1 : itermax 
        unew(I,I) = 0.25*(u(I,Ip) + u(I,Im) + u(Ip,I) + u(Im,I) - h2*f(I,I));
        ejac(k) = max(max(abs(unew - u_exact)));
        u = unew;
    end
    toc
    plot(ejac,'m','Linewidth',2);
    z = log(ejac(itermin : itermax));
    pol = polyfit(ii,z,1);
    fprintf('Jacobi, h = %d: convergence rate: %d\n',h(p - p1),-pol(1));
    rjac(p - p1) = -pol(1);

    % Gauss-Seidel
    u = zeros(n);
    tic
    for k = 1 : itermax
    for i = 2 : n1
        for j = 2 : n1
            u(i,j) = 0.25*(u(i,j+1) + u(i,j-1) + u(i+1,j) + u(i-1,j) - h2*f(i,j));
        end
    end
    egs(k) = max(max(abs(u - u_exact)));
    end
    toc
    plot(egs,'g','Linewidth',2);
    z = log(egs(itermin : itermax));
    pol = polyfit(ii,z,1);
    fprintf('Gauss-Seidel, h = %d: convergence rate: %d\n',h(p - p1),-pol(1));
    rgs(p - p1) = -pol(1); 
    
    % SOR
    u = zeros(n);
    om = 2/(1 + sin(pi*h(p - p1))) + eps;
    tic
    for k = 1 : itermax
    for i = 2 : n1
        for j = 2 : n1
            u(i,j) = (1 - om)*u(i,j) + om*0.25*(u(i,j+1) + u(i,j-1) + u(i+1,j) + u(i-1,j) - h2*f(i,j));
        end
    end
    esor(k) = max(max(abs(u - u_exact)));
    end
    toc
    plot(esor,'c','Linewidth',2);
    set(gca,'YScale','log') 
    z = log(esor(itermin : itermax));
    pol = polyfit(ii,z,1);
    fprintf('SOR, h = %d: convergence rate: %d\n',h(p - p1),-pol(1));
    rsor(p - p1) = -pol(1);
end
figure(2);
hold on;
plot(h,rjac,'r','Linewidth',2);
plot(h,rgs','b','Linewidth',2);
plot(h,rsor','k','Linewidth',2);
end


