P1 = @(x) x.^2 - x + 1/2;
P2 = @(x) -x.^2 + x + 1/2;
P2_P1 = @(x) P2(x) - P1(x);
sol = 1/3;

x = linspace(-0.5,1.5,1000);
figure(1)
plot(x,P1(x),'-')
hold on
plot(x,P2(x),'-')
hold off
legend('P1','P2')


N = 2.^(2:6);
M = length(N);
K = 500;
a = 0;
b = 1;
%% Type 1

error1 = zeros(M,1);

for i = 1:M
    error_N = zeros(K,1);
    for k = 1:K
        error_N(k) = abs(mc1(P2_P1,a,b,N(i)) - sol);
    end
    error1(i) = mean(error_N);
end
figure(2)
loglog(N,error1,'-')
hold on
loglog(N,1./sqrt(N),':')
hold off

%% Type 2

error2 = zeros(M,1);

for i = 1:M
    error_N = zeros(K,1);
    for k = 1:K
        error_N(k) = abs(mc2(P2_P1,a,b,N(i)) - sol);
    end
    error2(i) = mean(error_N);
end
figure(3)
loglog(N,error2,'-')
hold on
loglog(N,1./sqrt(N),':')
hold off

figure(4)
figure(3)
loglog(N,error1,'-')
hold on
loglog(N,error2,'-')
loglog(N,1./sqrt(N),':')
hold off
legend('MC1','MC2','1/sqrt(N)')
