S0 = 100;
sigma = 0.2;
rf = 0.05;
U = 120;
K = 105;

T = 1;

Ndt = 252; % number of steps from 0 to T
dt = T / Ndt;

M = 1000; % number of paths

S = NaN(M, Ndt + 1);
S(:,1) = S0;

Val =  zeros(M,1); % store the value of each scenario

for k = 1 : Ndt
   
    Z = randn(M,1);
    
    S(:, k+1)  = S(:, k) .* exp( (rf-0.5*sigma^2)*dt + sigma*sqrt(dt) * Z );
    
    idx = ( S(:,k+1) >= U );
    
    timeToMat = T-k*dt;
    dp = ( log(S(:, k+1) / K) + ( rf+0.5*sigma^2)*(timeToMat) ) / (sigma*sqrt(timeToMat));
    dm = dp - sigma*sqrt(timeToMat);
    

    % only store value the first time you hit the barrier
    idx2 = logical(idx.*(Val==0));
    Val(idx2) = ( S(idx2,k+1) .* normcdf( dp(idx2) ) - K*exp(-rf*timeToMat) * normcdf(dm(idx2)) ) *exp(-rf*k*dt); 
        
    
end

