function [P, CI] = MCAVButterfly(S0,r,T,sigma,NRepl,X1,X2,X3)
nuT = (r-0.5*sigma^2)*T;
siT = sigma*sqrt(T);
Veps = randn(NRepl,1);
Stocks1 = S0*exp(nuT + siT*Veps);
Stocks2 = S0*exp(nuT - siT*Veps);
Payoff1 = zeros(NRepl,1);
Payoff2 = zeros(NRepl,1);
In = find((Stocks1 > X1) & (Stocks1 < X2));
Payoff1(In) = (Stocks1(In) - X1);
In = find((Stocks1 >= X2) & (Stocks1 < X3));
Payoff1(In) = (X3 - Stocks1(In));
In = find((Stocks2 > X1) & (Stocks2 < X2));
Payoff2(In) = (Stocks2(In) - X1);
In = find((Stocks2 >= X2) & (Stocks2 < X3));
Payoff2(In) = (X3 - Stocks2(In));
Payoff = 0.5 * exp(-r*T) * (Payoff1 + Payoff2);
[P, V, CI] = normfit(Payoff);