function [u,loops] = PSOR_solver( u, b, g, a, omega, n_eps )
% PSOR_SOLVER - SOR algorithm for European options problems
%               this is a MATLAB vectorized version of the pseudo-code given 
%               in Figure 8.9 EPage 166 from the book The Mathematics of Financial Derivatives
%               by Paul Wilmott.
% 
% Input: 
%   u     = a vector of initial conditions 
%   b     = 
%   a     = deltaT/deltaX^2  
%   omega = a relaxation parameter 
%   n_eps = a convergence parameter 
%
% Output: 
%   u     = 
%   loops = 
% 
% Written by:
% -- 
% John L. Weatherwax                2008-06-11
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

loops = 0; 
error = 1e6; 
nelts = length(u); 

%fprintf('u(1) = %10.6f; g(1) = %10.6f\n', u(1), g(1) ); 
%fprintf('u(end) = %10.6f; g(end) = %10.6f\n', u(end), g(end) ); 
u(1)   = max( u(1), g(1) ); 
u(end) = max( u(end), g(end) ); 

while( error>n_eps )
  error = 0; 
 
  for n=2:(nelts-1),
    y     = ( b(n) + a * ( u(n-1) + u(n+1) ) )/(1+2*a); 
    y     = max( u(n) + omega * ( y - u(n) ), g(n) );
    error = error + (u(n)-y)^2;
    u(n)  = y;                                % <- use u(n) on the next iteration immediatly ... 
  end
    
  loops = loops+1; 

end