Как интегрировать в матлаб?

Q(aa)=integral2(@func1,0,2*AA+1,0,2*AA+1,'Method',mtd);

Q=integral2(@func1,0,2*AA+1,0,2*AA+1,'Method',mtd);

function z=zzz(kk,ll)
     global c a b z kk ll
     z=kk+ll
                  function kk=kkkk(c,a)
                          kk=c.*a;
                  end
                 function ll=llll(c,b)
                          ll=c./b;
                  end
end

global c a b z kk ll
c=1;
b=2;
a=3;

global A q Q QQQ AAA QQ AA
A=-2:1:2
q=-2:1:2

function F5zn = funcF5zn(A,q)
global A q Q QQQ AAA QQ AA
for aa=1:length(A)
    for bb=1:length(q)
        QQ=q(bb);
        AA=A(aa);
  QQQ=QQ.*cos(AA)+AA.*sin(AA);
 AAA=QQ.*sin(AA)+AA.*cos(AA); 
   zn=AA*cos(AA);
  F5zn (aa,bb) =dblquad(@funcB5n,0,2*AAA+1+QQQ+zn,0,2*AAA+1+QQQ);
    end
end
AAA
F5zn
function B5n = funcB5n(xn,bn)
  x=xn.*cos(AA)+bn.*sin(AA)+zn;
  b=bn.*sin(AA)+xn.*cos(AA);
      B5n = QQ.*AA.*(x./abs(x)).*(b./abs(b))+zn;
B5n( (b==0)|(x==0) ) = zn;
  end
end

AAA =
 
    0.9863
 
 
F5zn =
 
  415.9935  139.9977   19.8338   -4.4273    0.0560
    5.9904    0.0843   -0.0907   -0.3152   -3.3050
         0         0         0         0         0
   -0.9060   -0.3981    5.4664   45.0605  148.0412
   -8.9492    0.0438   -0.3090   -1.1657   16.3150

function F5zn2 = funcF5zn2(A,q)
global A q Q QQQ AAA QQ AA
    for aa=1:length(A)
    for bb=1:length(q)
      QQ=q(bb);
      AA=A(aa);
  QQQ(aa,bb)=QQ.*cos(AA)+AA.*sin(AA);
 AAA(aa,bb)=QQ.*sin(AA)+AA.*cos(AA); 
   zn=AA*cos(AA);
  F5zn2(aa,bb)=dblquad(@funcB5n2,0,2*AAA(aa,bb)+1+QQQ(aa,bb)+zn,0,2*AAA(aa,bb)+1+QQQ(aa,bb));
    end
end
AAA
 F5zn2
  function B5n2 = funcB5n2(xn,bn)
  x=xn.*cos(AA)+bn.*sin(AA)+zn;
  b=bn.*sin(AA)+xn.*cos(AA);
      B5n2 = QQ.*AA.*(x./abs(x)).*(b./abs(b))+zn;
B5n2( (b==0)|(x==0) ) = zn;
  end
end

AAA =
 
    2.6509    1.7416    0.8323   -0.0770   -0.9863
    1.1426    0.3012   -0.5403   -1.3818   -2.2232
         0         0         0         0         0
   -1.1426   -0.3012    0.5403    1.3818    2.2232
   -2.6509   -1.7416   -0.8323    0.0770    0.9863
 
 
F5zn2 =
 
  415.9935  139.9977   19.8338   -4.4273    0.0560
    5.9904    0.0843   -0.0907   -0.3152   -3.3050
         0         0         0         0         0
   -0.9060   -0.3981    5.4664   45.0605  148.0412
   -8.9492    0.0438   -0.3090   -1.1657   16.3150

clear global A q Q QQQ AAA QQ AA

global A q Q QQQ AAA QQ AA
A=-2:1:2
q=-2:1:2
F5zn2 = funcF5zn2(5,6)

function Q = powerF7yn(A)
Q = dblquad(@func1, 0,A,0,2)
 
function z = func1(x,y)
z = (x./abs(x))+(y./abs(y));
z((x==0))=0;
z((y==0))= 0;
end
end

function Q = powerF7yn(A)
Q = dblquad(@func1, 0,A,0,2)
 
function z = func1(x,y)
z = (x./abs(x))+(y./abs(y));
z((x==0)) = (y./abs(y));
z((y==0)) =(x./abs(x));
end
end

???  In an assignment  A(I) = B, the number of elements in B and
 I must be the same.
 
Error in ==> powerF7yn>func1 at 7
z((y==0)) =(x./abs(x));
 
Error in ==> quad at 77
y = f(x, varargin{:});
 
Error in ==> dblquad>innerintegral at 78
    Q(i) = quadf(intfcn, xmin, xmax, tol, trace, y(i), varargin{:});
 
Error in ==> quad at 77
y = f(x, varargin{:});
 
Error in ==> dblquad at 54
Q = quadf(@innerintegral, ymin, ymax, tol, trace, intfcn, ...
 
Error in ==> powerF7yn at 2
Q = dblquad(@func1, 0,A,0,2)

function Q = powerF7yn(A)
Q = dblquad(@func1, 0,A,0,2)
 
function z = func1(x,y)
z = (x./abs(x))+(y./abs(y));
z((x==0))=0;
z(repmat(y==0,size(x)))= 0;
end
end

function Q = powerF7yn(A)
Q = dblquad(@func1, 0,A,0,2)
 
function z = func1(x,y)
z = x./abs(x)+y./abs(y);
z(x==0)=0;
z(y==zeros(size(x)))= 0;
end
end

function Q = powerF7yn(A)
Q = dblquad(@func1, 0,A,0,2)
 
function z = func1(x,y)
z = (x./abs(x))+(y./abs(y));
z(x==0|y==0)=0;
end
end

function Q = powerF7yn(A)
Q = dblquad(@func1, 0,A,0,2)
 
function z = func1(x,y)
z = (x./abs(x))+(y./abs(y));
z((x==0)) = (y./abs(y)); % Вектору можно присвоить скаляр, 
%тогда все элемента вектора станут равны скаляру
z((y==0)) =(x./abs(x)); % Одному элементу нельзя присваивать вектор, 
% если этот елемент, конечно, не элемента массива ячеек
end
end

function Q = powerF9yn(A)
Q = 2*quadgk(@func1, 0,A)+A*quadgk(@func2, 0,2)
 
function z = func1(x)
z = x./abs(x);
% z(x==0) = 0;
end
 
function z = func2(y)
z = y./abs(y);
% z(y==0) =0;
end
 
end

function Q = powerF7yn(A,a)
Q = dblquad(@func1, 0,A,0,2)
 
function z = func1(x,y)
z = (x./abs(x))+(y./abs(y))+a./abs(a);
z((x==0))=0;
z(repmat(y==0,size(x)))= 0;
z(repmat(a==0,size(x)))= 0;
end
end

Q=trapz(x,y);

x=[0 1];
y=[1 1];

function Q = powerF15yn(x,y)
% pp=pchip(x,y); %Можно и полиномами Эрмита
pp=spline(x,y); %Можно и сплайном
Q=diff(fnval(fnint(pp),[min(x), max(x)]));
end

function Q = powerF16yn(x,y)
% aFittype = fittype('a*x+b')
aFittype = fittype('exp1') %как-то определился какой функцией апроксимировать
fitobject = fit(x(:),y(:),aFittype) %выполнил апроксимацию
Q = integral(@func1, min(x),max(x));
  function Result=func1(x)
  Result = feval(fitobject,x).';
  end
end

function Q = powerF(D,an,A,B,Y)
Q = dblquad(@func1, 0,D,0,2)
 
function z = func1(xn,yn)
x=xn.*cos(A).*cos(B) - yn.*sin(A).*cos(B) + an.*sin(B);
y= xn.*(cos(A).*sin(B).*sin(Y)+sin(A).*cos(Y)) + yn.*( -sin (A).*sin(B).*sin(Y)+ cos(A).*cos(Y) ) - an.*cos(B).*sin(Y);
a= -xn.*( cos(A).*sin(B).*cos(Y)- sin(A).*sin(Y)) + yn.*( sin (A).*sin(B).*cos(Y) + cos(A).*sin(Y) ) + an.*cos(B).*cos(Y);
 
z = (x./abs(x))+(y./abs(y))+a./abs(a);
z((x==0))=0;
z(repmat(y==0,size(xn)))= 0;
z(repmat(a==0,size(xn)))= 0;
end
end

function Q = powerF(D,an,A,B,Y)
Q = dblquad(@func1, 0,D,0,2)
 
function z = func1(xn,yn)
yyn=repmat(yn,size(xn));
aan=repmat(an,size(xn));
 
x=xn.*cos(A).*cos(B) - yyn.*sin(A).*cos(B) + aan.*sin(B);
y= xn.*(cos(A).*sin(B).*sin(Y)+sin(A).*cos(Y)) + yyn.*( -sin (A).*sin(B).*sin(Y)+ cos(A).*cos(Y) ) - aan.*cos(B).*sin(Y);
a= -xn.*( cos(A).*sin(B).*cos(Y)- sin(A).*sin(Y)) + yyn.*( sin (A).*sin(B).*cos(Y) + cos(A).*sin(Y) ) + aan.*cos(B).*cos(Y);
z = (x./abs(x))+(y./abs(y))+a./abs(a);
z((x==0))=0;
z((y==0))= 0;
z((a==0))= 0;
 
end
end

function Q = powerF1(D,an,A,B,Y)
Q = dblquad(@func1, 0,D,0,2)
 
function z = func1(xn,yn)
yyn=repmat(yn,size(xn));
aan=repmat(an,size(xn));
 
x=xn.*cos(A).*cos(B) - yyn.*sin(A).*cos(B) + aan.*sin(B);
y= xn.*(cos(A).*sin(B).*sin(Y)+sin(A).*cos(Y)) + yyn.*( -sin (A).*sin(B).*sin(Y)+ cos(A).*cos(Y) ) - aan.*cos(B).*sin(Y);
a= -xn.*( cos(A).*sin(B).*cos(Y)- sin(A).*sin(Y)) + yyn.*( sin (A).*sin(B).*cos(Y) + cos(A).*sin(Y) ) + aan.*cos(B).*cos(Y);
z = (x./abs(x))+(y./abs(y))+a./abs(a);
z((x==0))=0;
z((y==0))= 0;
z((a==0))= 0;
 
end
end
 
Q=0
Q=12

function Q = powerF2(D,a)
Q = dblquad(@func1, 0,D,0,2)
 
function z = func1(x,y)
 
 z = (x./abs(x))+(y./abs(y))+a./abs(a);
z((x==0))=0;
z(repmat(y==0, size(x)))= 0;
z(repmat(a==0, size(x)))= 0;
 
end
end
 
Q=0
Q=12

x=xn.*cos(A).*cos(B) - yn.*sin(A).*cos(B) + an.*sin(B);
y= xn.*(cos(A).*sin(B).*sin(Y)+sin(A).*cos(Y)) + yn.*( -sin (A).*sin(B).*sin(Y)+ cos(A).*cos(Y) ) - an.*cos(B).*sin(Y);
a= -xn.*( cos(A).*sin(B).*cos(Y)- sin(A).*sin(Y)) + yn.*( sin (A).*sin(B).*cos(Y) + cos(A).*sin(Y) ) + an.*cos(B).*cos(Y);

yyn=repmat(yn,size(xn));
aan=repmat(an,size(xn));

clear
x=0:0.1:10;
for xx=1:length(x)
y(xx)=x(xx).^2;
end
 
for zz=2:length(x)
z(zz)=trapz(0:0.1:x(zz),[B][I][U]y(zz)[/U][/I][/B])
end

clear
x=0:0.1:10;
for xx=1:length(x)
  y(xx)=x(xx).^2;
end
 
for zz=2:length(x)
  z(zz)=trapz(x(1:zz),y(1:zz))
end

x = 0:0.1:10;
y = x.^2;
z= cumtrapz(x,y);
plot(x,z)

clear
F=[0 0 0; 1 1 1; 2 4 8; 3 9 27; 4 16 64];
p=0:4;
for ff=2:length(p)
  A(ff)=trapz(p(1:ff), F(1: ff));
end
figure 
plot(p, F )
figure
plot(p, A )

clear
F=[0 0 0; 1 1 1; 2 4 8; 3 9 27; 4 16 64];
p=0:4;
A = cumtrapz(p, F );
 
figure 
plot(p, F )
figure
plot(p, A )