Integrand output size does not match the input size

Hello Guys!!! I am working in this code for hours, I undertand the rule of the integral2

"The function fun must accept two arrays of the same size and return an array of corresponding values. It must perform element-wise operations." and usually when I get a error like that is some ./ missing... but for me know everything is correct, I dont understand what is going on in this code

%dimensions

d1=40;

d2=d1;

d=d2;

y=sind(45).*d;

z=y;

alfa=0;

%triangle 1 quadrants 2 and 1

dy1=0;

dz12=0;

dz11=0;

%triangle 2 quadrants 2 and 3

dz2=0;

dy22=0;

dz23=0;

%triangle 3 quadrants 3 and 4

dy3=0;

dz33=0;

dz34=0;

%triangle 4 quadrants 4 and 1

dz4=0;

dy44=0;

dz41=0;

%%%%%%%%%%area quadrant 2

fun = @(y,z) 1;

%Area 2B

y2bmax = @(z) sqrt(d + (z.^2));

y2bmin = 0;

z2bmax = z - dz2;

z2bmin = dz12;

a2b = integral2(fun,z2bmin,z2bmax,y2bmin,y2bmax);

z2b = integral2(z,z2bmin,z2bmax,y2bmin,y2bmax);

y2b = integral2(y,z2bmin,z2bmax,y2bmin,y2bmax);

%Area 2A

y2amax = y-dy1;

y2amin = 0;

z2amax = dz12;

z2amin = 0;

a2a = integral2(fun,z2amin,z2amax,y2amin,y2amax);

z2a = integral2(z,z2amin,z2amax,y2amin,y2amax);

y2a = integral2(y,z2amin,z2amax,y2amin,y2amax);

at2=a2b+a2a;

z2=(z2b+z2a)./(at2);

y2=(y2b+y2a)./(at2);

f2=1.8.*(4.56.*(10.^(-6))).*((cos(2*alfa)+1)./2)*(at2);

Error using integral2Calc>integral2t/tensor (line 242)

Integrand output size does not match the input size.

Error in integral2Calc>integral2t (line 56)

[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

Error in integral2Calc (line 10)

    [q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral2 (line 107)

    Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);

Error in my_ode>EOM (line 56)

a2b = integral2(fun,z2bmin,z2bmax,y2bmin,y2bmax);

Error in odearguments (line 88)

f0 = feval(ode,t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode45 (line 114)

[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in my_ode (line 5)

    [t,o]=ode45(@EOM,[0 100],[0,0,0]);