Evaluating a triple integral in Matlab using the trapezoidal rule

Question

Can someone please explain to me how evaluate a triple integral in matlab using the trapezoidal rule.

I have to evaluate a triple integral with six digits of accuracy and these limits:

$0\leq x \leq 2 $

$0\leq y \leq sin(pi*x/2) $

$0\leq z \leq 1 $

$xh=1.25;$

$yh=0;$

$zh=0.25;$

This is the function:

$ f = exp(-(50*(z+zh)*exp(-z/2*zh))*(x-xh).^2+y.^2)*cos(y*(x-xh))$

Is my code right? When I run it it runs forever...

My Matlab Code:

function f=funcz(x,y,z)

    xh=1.25;

    yh=0;

    zh=0.25;

    f = exp(-(50*(z+zh)*exp(-z/2*zh))*(x-xh).^2+y.^2)*cos(y*(x-xh));

end


    clear all, close all, clc
    format long 

    %Constants given by the problem
    xh=1.25;
    yh=0;
    zh=0.25;

    %Trapezoidal rule
    %Step length
    hx=0.05;
    hy=0.05;
    hz=0.05;

    %Empty vectors 
    yt1 = [];
    Itraps1 = [];
    Itraps2 = [];
    Istep = [];
    error = [];    
    hxx = [];
    hyy = [];
    Er_step = 1;

    %Iteration counter
    i=1;
    k=1;
    j=1;

    %%

    %Matlab solution
    intergrand = @(x,y,z) exp(-(50.*(z+zh).*exp(-z/2.*zh)).*(x-xh).^2+y.^2).*cos(y.*(x-xh));
    y_max = @(x) sin(pi*x/2);
    Irekt_matlab = integral3(intergrand,0,2,0,y_max,0,1)

    %%

    while Er_step > 1.e-1
        xx=0:hx:2;    
        zz=0:hz:1;
        for z = 0:hz:1
            for x = 0:hx:2

                % if clause to avoid hy=0
                if x == 0 
                    hy = hy;
                elseif x == 2
                    hy = hy;
                else
                    yy=sin(pi*x/2);
                    hy = yy/round(yy/hx);
                end
                % if clause end

                for y = 0:hy:sin(pi*x/2)
                    f = funcz(xx(i),y,zz(j));
                    yt1 = [yt1 f]; 
                end
                Itraps = hy*(sum(yt1(2:end-1))+1/2*yt1(1)+1/2*yt1(end));
                Itraps1 = [Itraps1 Itraps];
                yt1 = []; 
                i=i+1;
            end
            Itraps = hx*(sum(Itraps1(2:end-1))+1/2*Itraps1(1)+1/2*Itraps1(end));
            Itraps2 = [Itraps2 Itraps];
            j=j+1;
            i=1;
        end
        I = hz*(sum(Itraps2(2:end-1))+1/2*Itraps2(1)+1/2*Itraps2(end));

        %Save integrals value after each step
        Istep = [Istep I];
        Itraps1 = [];
        Itraps2 = [];
        %Calculate the error after each step and save the value
        Er_step = abs(Irekt_matlab - Istep(k));
        error = [error Er_step];
        %Half the step and save the value
        hxx = [hxx hx];
        hyy = [hyy hy];
        hx = hx/2;
        hz = hz/2;
        hy = hx;
        i=1;
        j=1;
        k=k+1;
    end

    Isolution = Istep(end)

    p1 = log2((Istep(end-1)- Irekt_matlab)/(Istep(end)- Irekt_matlab))
    p2 = log2((Istep(end-2)- Istep(end-1))/(Istep(end-1)- Istep(end)))

    loglog(hxx,hxx.^2,'b')
    hold on
    loglog(hxx,error,'r')

    grid on