function [H,r,t] = IAR3D(alpha,beta,vol,expry,fixedrate,period,NRS,NTS)
% IAR3D - solve the index amortizing rate swap equation.
%   
% Written by:
% -- 
% John L. Weatherwax                2006-12-31
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

dt = expry/NTS;
dr = 3*fixedrate/NRS;

% compute the grid of interest rates: 
r = (0:NRS)' * dr;
% compute the grid of timesteps: 
t = (0:NTS)' * dt; 

% Compute the grid of final conditions on H:
%
% We integrate by marching along each column from the 
% first column (final condition) to the last (initial condition)
%
% H(:,1)   = is the final time t=T for all r ... 
% H(:,end) = is the initial time t=0 for all r 
% H(1,:)   = is the smallest interest rate r=0
% H(end,:) = is the largest interest rate 
H = zeros(NRS+1,NTS+1);
H(:,1) = r - fixedrate;

for k=1:NTS, % column index 1=>t=T, end=>t=0 
  tim = (k-1)*dt; % = time till expiration = T-t 
  delta = ( H(3:(NRS+1),k) - H(1:(NRS-1),k) ) / (2.0*dr); 
  gamma = ( H(3:(NRS+1),k) - 2 *  H(2:NRS,k) + H(1:(NRS-1),k) )/( dr*dr );
  theta = -0.5 * vol * vol * r(2:NRS) .* r(2:NRS) .* gamma - ...
          ( alpha - beta * r(2:NRS) ) .* delta + r(2:NRS) .* H(2:NRS,k); % bond pricing equation 
  H(2:NRS,k+1) = H(2:NRS,k) - dt * theta; 
  % boundary contitions: 
  slopey       = (H(2,k) - H(1,k))/dr;
  timederiv    = -alpha * slopey;
  H(1,k+1)     = H(1,k) - dt * timederiv;
  H(NRS+1,k+1) = 2 * H(NRS,k+1) - H(NRS-1,k+1);
  % jump condition: 
  if( floor(tim/period) < floor((tim+dt)/period) - 0.001 )
    H(:,k+1) = amortizing_schedule(r) .* H(:,k+1) + r - fixedrate; 
  end
end