fprintf('\n\n\n\n\n\n\n\n\n\n\n\n\n');
points=input('Please enter number of points in the 1m gap '); % Get user input for amount of points
iterations=input('Please enter maximum number of iterations '); % Get user input for max iterations

Nx=points+1; Ny=Nx; h=1/Nx; %Set h and Nx
thi=zeros(points+2,points+2);R=thi; % Setup matrices with all zeros
thi(1,1)=5; thi(1,2:Nx)=20; thi(1,Nx+1)=15; thi(2:Ny,1)=-10; thi(Ny+1,1)=-5; thi(2:Ny,Nx+1)=10;
thi(Ny+1,Nx+1)=5; thi(Ny+1,2:Nx)=0; thi(2:Ny,2:Nx)=5;% Setup initial values of thi
thinew=thi; % Have thinew the same as thi initialised

fprintf('\nStats:\n')
fprintf('distance between points %f m\n',h)

a=10^-9; E=8.8542*10^-12; % Set constants
t=cos(pi/Nx) + cos(pi/Ny); % Evaluate t for matrix
temp=[t^2 -16 16]; ws=roots(temp); %Find the roots of equation to find values for w
if ws(1)>ws(2)
    w=ws(2);
else
    w=ws(1);
end

tic; % Measure time taken
for loop=1:iterations
    error=0; %Set error to 0 at beginning of each iteration
    for i=2:Nx
        for j=2:Ny
            x=(i-1)*h; %Find distance in x
            y=1-(j-1)*(h); %Find distance in y
            g(i,j)=-a*x*(y-1)/E; % Evaluate charge equation
            R(i,j) = thi(i+1,j) + thi(i-1,j) + thi(i,j+1) + thi(i,j-1) - 4*thi(i,j) - (h^2)*g(i,j); % Find Residual error matrix
            thinew(i,j)=thi(i,j) + (w/4)*R(i,j); % Ierative step to update thi with respect to residual
            error=error+abs( thinew(i,j) - thi(i,j) ); %Evaluate an error for stopping if converged
            thi=thinew;
        end
    end
    
    if error<(10^-5)/points^2 % If small enough error then converged
        break % Therefore break the iterations
    end
    
end

%Display results
timetaken=toc;
fprintf('After %d iterations and %f seconds, the final result is:',loop,timetaken)
thi(2:Nx,2:Ny)

%Plot data in 3D graph
x=linspace(0,1,Nx+1);
y=x;
[X,Y]=meshgrid(x,y);
surf(X,Y,thi);
xlabel('x (metres)')
ylabel('y (metres)')
zlabel('charge (coulombs)')
title('Charge distribution')