function u = lu_solver( u, b, y, a )
% See the Epage 163 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.
% 
%-----

q = zeros(size(u)); 

Nminus = 1; 
Nplus  = length(y); 

q(2) = b(2);
%q(3:(end-1)) = b(3:(end-1)) + a*q(2:(end-2))./y(2:(end-1)); % <- incorrect .. can't vectorize in this way 
for n=3:(Nplus-1), 
  q(n) = b(n) + a*q(n-1)/y(n-1);
end

u(Nplus-1) = q(Nplus-1)/y(Nplus-1); 

for n=(Nplus-2):-1:(Nminus+1), 
  u(n) = (q(n)+a*u(n+1))/y(n); 
end