function [u,xgrid] = crank_fd_SOR( payoff_fn, u_m_inf, u_p_inf, r, sigma, xLeft, xRight, Nx, tau_Max, M )
% 
% Solves the diffusion equation using the Crank-Nicolson finite-difference model using the SOR solver
% 
% See the Epage 177 in Wilmotts book: The Mathematics of Financial Derivatives
% 
% 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.
% 
%-----

xgrid = linspace( xLeft, xRight, Nx+1 ); 
dx    = (xRight-xLeft)/Nx; 
dt    = tau_Max/M; 

a        = dt/(dx*dx); 
a2       = 0.5*a; 
n_eps    = 1.e-8; 
omega    = 1.0; 
domega   = 0.05; 
oldloops = 10000; 
k        = r/(0.5*sigma^2); 

% get the initial conditions: 
tau  = 0.0; 
oldu = payoff_fn( xgrid, tau, k ); oldu = oldu(:).'; 

u      = zeros(M+1,Nx+1);    % <- space for our results ...
u(1,:) = oldu; 
newu   = zeros(1,Nx+1); 
b      = zeros(1,Nx+1); 
for m=1:M,
  tau = m*dt; 

  b( 2:(end-1) ) = (1-a)*oldu( 2:(end-1) ) + a2*( oldu( 3:end ) + oldu( 1:(end-2) ) );
  
  % update the endpoints: 
  newu(1)   = u_m_inf( xgrid(1),   tau, k ); 
  newu(end) = u_p_inf( xgrid(end), tau, k ); 

  [newu,loops] = SOR_solver( newu, b, a2, omega, n_eps );
  if( loops > oldloops ) domega = -domega; end; 
  omega    = omega + domega; 
  oldloops = loops; 
  % prepare for the next pass: 
  oldu = newu;   

  % store 
  u(m+1,:) = newu; 
end