Svp Aide code Matlab

Bonjour,

j'ai un code Matlab pour 1 dimension, j'ai besoin de votre aide pour écrit même code pour 2 dimension.

J’ai besoin de montrer que l'erreur Monte-Carlo est fonction de 1/sqrt(N) . Mon prof me dit que : Pour obtenir les erreurs MC, pour chaque nombre de points N (10^2, 10^3, 10^4...10^7) vous refaites le calcul (30 fois par exemple) pour calculer sigma, et l'erreur sera sigma/sqrt(N). vous aurez l'erreur pour chaque N, vous tracer la courbe pour obtenir la courbe en bleu. Et tracer 1/sqrt(N) qui est la courbe verte. Il me donne un exemple d’intégrale pour 1D ( code+commentaires explicatif ) et il me demande de faire la même chose pour 2D, pour montré que l’erreur converge vers 1/sqrt(N) n'importent qu’il dimension. Mais je n’arrive pas affaire , j’ai besoin l'aide svvvp.

code 1D :

close all
clc
format long

N = [10^2 10^3 10^4 10^5 10^6 10^7]';
C = 0.25; %fmax
a = 0;
b = 1;

%Nombre des essais
MCTries = 7;

%On aura besoin de ces 2 tableaux pour garder enregister lerreur 
%correspondante à chaque essaye (pour cet exemp: pour chaque N on refera le calcul 7 fois.
%on prendera à la fin la moyenne des erreurs.
array_sigmas = zeros(MCTries,1);
array_errors = zeros(length(N),1);

for i=1:length(N) % boucle sur les N
    for j=1:MCTries % boucle sur les essais
        I = 0;
        f = zeros(N(i,1),1);
        for k = 1:N(i,1) % boucle sur les xi
            x = rand();
            fx = x*(1-x)*sin(200*x*(1-x))^2;
            I = I + fx;
            f(i,1) = fx;%stockage des fx
        end
        % I = (b-a)*I/N(i,1) % valeur de l'integrale
    array_sigmas(j,1) = sqrt(var(f)); % calcul de l'erreur (la racine de la variance)  
    end
    % on prends la moyenne des essais
    array_errors(i,1) = mean(array_sigmas);
end


%data to plot
 x = N(:,1);
y0 = 1./sqrt(N(:,1));
y1 = array_errors(:,1)

%plotting
figure1 = figure;
axes1 = axes('Parent',figure1);
hold(axes1,'on');
box(axes1,'on');
grid(axes1,'on');
plot(x,y0,'-r',x,y1,'-bo');
set(axes1,'XMinorTick','on','XScale','log','YMinorTick','on','YScale','log');
xlabel('N') % x-axis label
ylabel('errors') % y-axis label
legend('1/sqrt(N)','error in MC integration')

pour 2D : par exemple cette fonction f(x,y)=xye^((-x^2)y) ou bien une fonction de votre choix.
Et merci