clear all;close all

Nx=20;Ny=20;                             % Number of points
Lx=1;Ly=1;                               % Length of the domain  
omega=1.0;                               % SOR parameter
hx=Lx/Nx;hy=Ly/Ny;                       % Grid spacing
beta = hx/hy; b2 = beta^2; hx2=hx^2;     % Grid spacing parameter   
x=linspace(0,Lx,Nx+1);                   % x-coordinates length N+1 at grid nodes
y=linspace(0,Ly,Ny+1);                   % y-coordinates length N+1 at grid nodes

% Assemble x and y - coordinate at cell centers into 1-D vector length N*N
X=(x(1:end-1)+x(2:end))/2;
Y=(y(1:end-1)+y(2:end))/2;

% Define the forcing function (NOTE: correct ordering!)
f=...

%-----------------------------------------
Lp = kron(speye(Ny),...) + kron(...,speye(Nx));
% Set one Dirichlet value to fix pressure in that point
Lp(1,:) = ...; Lp(1,1) = ...;f(1)=...;

P = LP\f; 

figure(2)
surf(X,Y,reshape(P,Nx,Ny)')
zlabel('P_{dir}');xlabel('x');ylabel('y');title('Direct Solution')
%-----------------------------------------

% Set up u for iterative method
p=zeros(Nx*Ny,1)    ;                         % Inital condition

% Choose a criterion for convergence and the maximum iteration steps
err=1e-12;
Mmax=10000;
m=0;
convergence=1;

while convergence>err &&  m<Mmax
    m=m+1;
    
    % B.C. P_[1,1]
    p(1)=...
    
    for i=2:(Nx-1)
        % B.C. P_[i,1]
        p(i)=...
    end
    
    % B.C. P_[Nx,1]
    p(Nx)=...
    
    for j=2:(Ny-1)
        % B.C. P_[1,j]
        k=1+Nx*(j-1);
        p(k)=...
        % General Solution (March through all inner points)
        for i=2:(Nx-1)
            k=i+Nx*(j-1);
            p(k)=...
        end
        % B.C. P_[Nx,j]
        k=Nx+Nx*(j-1);
        p(k)=...
    end
    
    % B.C. P_[1,Ny]
    k=1+Nx*(Ny-1);
    p(k)=...
    
    for i=2:(Nx-1)
        % B.C. P_[i,Ny]
        k=i+Nx*(Ny-1);
        p(k)=...
    end
    
    % B.C. P_[Nx,Ny]
    k=Nx*Ny;
    p(k)=...
    
    
    % Calculate residuals
    R(m)=...
    convergence=...
    
    % Calculate difference to exact solution
    epsilon(m)=...
    
end

figure(3)
surf(X,Y,reshape(p,Nx,Ny)');
xlabel('x');ylabel('y');zlabel('P_{iter}');zlim([0 0.1]);
title(['Iterative solution after steps=',num2str(m)])

figure(4);
semilogy(R);hold on
xlabel('Iterations')
ylabel('log_{10}(R)')

figure(5);
semilogy(epsilon);hold on
xlabel('Iterations')
ylabel('log_{10}(\epsilon)')

% slope=(((log10(R(9000))-(log10(R(12000)))))/(12000-9000))

