function [P, CI] = MCButterfly(S0,r,T,sigma,NRepl,X1,X2,X3)
nuT = (r-0.5*sigma^2)*T;
siT = sigma*sqrt(T);
Veps = randn(NRepl,1);
Stocks = S0*exp(nuT + siT*Veps);
In1 = find((Stocks > X1) & (Stocks < X2));
In2 = find((Stocks >= X2) & (Stocks < X3));
Payoff = exp(-r*T)*[(Stocks(In1)-X1);  (X3-Stocks(In2)); ...
         zeros(NRepl - length(In1) - length(In2),1)];
[P, V, CI] = normfit(Payoff);