function [ S AV alpha beta time Boundary ] = American(S0, K, T, sigma, rf, Nsteps)

    % compute the various model parameters
    
    % time step size
    dt = T/Nsteps;
    
    % risk-neutral probabilities
    u = exp(sigma*sqrt(dt));
    q = ( exp(rf*dt) - 1/u ) / (u-1/u);
    
    disc = exp(-rf*dt);

    %
    % initialize payoffs, strategy and exercise boundary
    %
    AV = NaN(Nsteps+1, Nsteps+1);
    S = NaN(Nsteps+1,Nsteps+1);
    alpha = NaN(Nsteps+1, Nsteps+1);
    beta = NaN(Nsteps+1, Nsteps+1);
    Boundary = NaN * ones(Nsteps+1,1);
    time = zeros(Nsteps+1,1);
    
    %
    % set payoff and boundary at maturity    
    %
    S(:,end) = S0 * u.^(Nsteps - 2*[0:Nsteps])';
    AV(:,end) = (S(:,end) < K) .* (K-S(:,end));

    Boundary(Nsteps+1) = K;
    time(Nsteps+1)=T;
    
    %
    % step backwards through time
    %
    % "i" records number of steps *from* maturity
    % so that i = 1 refers to step n-1
    %
    for i = 1 : Nsteps
        
        n = Nsteps+1-i;
        %
        % step down through the tree starting from the top most node at the
        % current time step
        %
        AV([1 : n ], n) = disc*( q * AV([1:n], n+1) + (1-q)* AV([2:n+1], n+1) );
        
        % update spot prices...
        S([1:n],n) = S([1:n],n+1)/u;

        % find optimal exercise point...
        idx = find( AV([1:n],n) < K-S([1:n],n), 1, 'first');
        if ~isempty(idx)
            Boundary(n) = S(idx,n);
        end        
        AV([1:n],n) = max( AV([1:n],n), K-S([1:n],n));
        
        % find replicating portfoli...
        alpha([1:n], n) =     ( AV([1:n], n+1) - AV([2:n+1], n+1)) ./ (S([1:n],n+1) - S([2:n+1],n+1));
        beta([1:n], n) =  (1/disc^(n-1))*(AV([1:n], n) - S([1:n],n).*alpha([1:n],n));
        
        
        % record current time step
        time(Nsteps-i+1) = dt*(Nsteps-i);
        
    end
   
    
end