% Parameters
T = 1; S0 = 10; K = 10;
r = 0.01; sigma = 0.5; mu = 0.04;
Payoff = @(S) max(S-K,0);

% Barrier (down and out call)
B = 6;
% Black-Scholes value
S = linspace(B,21,201);
V = @(S) bs(S,K,r,sigma,T,0,'call') - ...
    ((S/B).^(1-2*r/sigma^2)).*bs(B^2./S,K,r,sigma,T,0,'call');
plot(S,V(S));

M = 20;
myerror = []; dist = [];
for j = 1:12
    j
    M = 2^j;
    [u,d,p,discount] = maketree(T,sigma,r,mu,M,1);
    % Final asset values
    S = S0*(u.^(0:1:M)).*(d.^(M:-1:0));
    [dist(j),i] = min(abs(S-B));
    BB = B;
    VB = Payoff(S);
    VB(S<=BB) = 0; % set values below the barrier to zero
    % loop back through time
    for m = M:-1:1
        S = S0*(u.^(0:1:m)).*(d.^(m:-1:0));
        VB = discount*(p*VB(2:end) + (1-p)*VB(1:end-1));
        VB(S<=BB) = 0; % set values below the barrier to zero
        %disp(VB)
    end
    disp('The binomial value is')
    disp(VB)
    disp('The Black-Scholes value was')
    disp(V(S0))
    myerror(j) = VB-V(S0);
end