%% Poisson Test Program
   clear all %clear all variables
   clc % clear output screen
   
   nx = 21; % Matlab's index for matrices starts from 1
   ny = nx; %different nx & ny gives non-square  
   w	= zeros(nx,ny); % form a matrix of zeros (1..nx, 1..ny) 
   p = w; % w = voriticty, p = streamfunction

   h = 1/(nx-1); %h = dx (=dy),
   r = 2/(1+pi/(nx-1)); % r = SOR parameter
   %so x-size is 1, y-size is ny/nx
   
   %% Initiaze the values for p and w
   for i = 1:nx
      for j =1:ny
	      p(i,j) = 0;
	      w(i,j) = 2*pi*pi*sin(pi*(i-1)*h)*sin(pi*(j-1)*h);
      end
   end

   %%Iterations start here
   
   for it = 1:4*nx 

      er2 = 0;
      for i = 2:nx-1 
	        for j = 2:ny-1 
	          q = p(i-1,j) + p(i+1,j) + p(i,j-1) + p(i,j+1)...
	              -4*p(i,j) + h*h*w(i,j) ;
	          p(i,j) =  p(i,j) +  r*0.25*q;
	          if (er2<abs(q)) 
              er2 = abs(q); %er2 = residue in satisfying equation
            end
          end
      end

      er = 0;
      for i = 1:nx
	      for j = 1:ny
	        if (abs(p(i,j)-sin(pi*(i-1)*h)*sin(pi*(j-1)*h)) > er)
	          er = abs(p(i,j)-sin(pi*(i-1)*h)*sin(pi*(j-1)*h));
            %er = error cf analytic answer
          end
        end
      end
      %Output of iterations
      formatSpec1='it=%d; er=%12.7f; er2=%12.7f; p() =%10.5f\n ';
      fprintf(formatSpec1, it, er,er2, p(floor(nx/2),floor(nx/2)));

   end
   
   
   %%Plot (3D) solution
   
   
   tx = linespace(0, (nx-1)*h, nx);
   ty = linespace(0,(nx-1)*h,nx);
   [xx, yy] = meshgrid (tx, ty);
   mesh (tx, ty, p);
   xlabel ('x');
   ylabel ('y');
   zlabel ('phi(x,y)');
   title ('Poisson Plot');
   



   