clear all 
% computational domain 
xmin = 0; xmax = 1. ;
% Discretization  and data 
N = 41; 
Big = 3 ;
x = linspace(xmin, xmax, N);
y = x;
h = x(2)-x(1);
[X,Y] = meshgrid(x,y) ; 
n = refract(X,Y) ; 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% init 
phi = Big *ones(N,N) ; 

% BC Dirichlet 
phi(1:N, 1) =  Big  ; 
phi(1:N, N) = Big ; 
phi(1, 1:N) = Big ; 
phi(N, 1:N)= Big ; 
%%% 
 
[phi,exact] = source(X,Y,phi,0.25,0.25,0.75,0.75,h) ; %% enforce source point 


%% iteration parameters 
conv = Big ;
tol = 1e-14 ; 
niter = 0;

phiold = phi; %% sequence 

%% firsr order ..

while ( conv > tol)  %% iterations 
niter = niter + 1
%%  

phi = sweep(phiold,n,h,N,[2:N-1],[2:N-1]) ; 

phiold = phi ; 

phi = sweep(phiold,n,h,N,[2:N-1],[N-1:-1:2]) ; 

phiold = phi ; 

phi = sweep(phiold,n,h,N,[N-1:-1:2],[2:N-1]) ; 

phiold = phi ; 

phi = sweep(phiold,n,h,N,[N-1:-1:2],[N-1:-1:2]) ; 

phiold = phi ; 

conv  = max(max(abs(phiold-phi)))/max(max(abs(phiold))) 
err  = max(max(abs(exact(2:N-1,2:N-1)-phi(2:N-1,2:N-1))))/max(max(abs(exact(2:N-1,2:N-1)))) 

phiold = phi  ; 

contour(x(2:N-1), y(2:N-1), phi(2:N-1,2:N-1)) 
colorbar
pause 

end 

contour(x(2:N-1), y(2:N-1), phi(2:N-1,2:N-1), 20,'b-') 


pause 
contour(x(2:N-1), y(2:N-1), exact(2:N-1,2:N-1), 20,'r-') 
pause 
contour(x(2:N-1), y(2:N-1), abs(phi(2:N-1,2:N-1)-exact(2:N-1,2:N-1))) 
colorbar