clear

%define a 21x21 matrix 
NUM=21;
THRESH=1e-6;
M=zeros(NUM);

%define a matrix which sums the nearest neighbors along one dimension and divides by 4
AVG=(diag(ones(NUM+1,1),1)+diag(ones(NUM+1,1),-1))/4;

newM=rand(NUM); %initialize
while (sum(sum(abs(newM-M))) > THRESH) %do the following until it doesn't make any difference to the outcome
     newM=M;%duplicate matrix

     %surround matrix with zeros at boundaries
     M=[zeros(1,NUM); M; zeros(1,NUM)];
     M=[zeros(NUM+2,1) M zeros(NUM+2,1)];
     
     M(12,12)=1; %set boundary condition at center point always equal to one: a point charge!
     
     M=AVG*M+(AVG*M')'; sum along x- and y-directions, dividing by 4

     M=M(2:(NUM+1),2:(NUM+1)); %discard zeros
end

%add boundaries again
M=[zeros(1,NUM); M; zeros(1,NUM)];
M=[zeros(NUM+2,1) M zeros(NUM+2,1)];
M(12,12)=1;

%plot potential
surf(M)
