Whats wrong if this triple integral ? Many errors output

I have a double integral,

      %seção de choque difusa

function exem6
    tic
f=@(b) anggbsR(b);
%x0=[1,2];
 x=fzero(f,1.1);
T=100;
g=10000;

f1=@(b) 2*(1 + (sqrt(pi*T)/(2*g)).*sin(ang2(b)/2)).*b ;
f2=@(b)  2*(1 - cos(ang2(b))).*b ;
T=100;

%qQd = integral (f2,x, 10,'RelTol',0,'AbsTol',1e-7) + integral(f1,0,x) 
qQd =  + integral(f1,0,x) + integral (f2,x, 10)
toc

      end

      function xx2 = ang2(b)
      for i=1:length(b)
          xx2(i)=anggbsR(b(i));
      end
      end

function xx=anggbsR(b)
g=10000;
T=100;
s= 0.6; 
R= 1; 
syms r
func =  @(r)  1 - (b/r)^2 - (g^-2)*((2/15)*[(s/R)]^9 *(1/(r - 1)^9 - 1/(r + 1)^9 - 9/(8* r) *(1/(r - 1)^8 - 1/(r + 1)^8)) - [(s/R)]^3 *(1/(r -1)^3 - 1/(r + 1)^3 - 3/(2* r)* (1/(r - 1)^2 - 1/(r + 1)^2)));
minrootgbsR = fzero(func,[1.01,100]);
k=1;
for i=1:length(minrootgbsR)
    if(isreal(minrootgbsR(i))==1)
        vec(k)=minrootgbsR(i);
        k=k+1;
    end

end

fun =@(r) 1./(r.^2 .* sqrt(1 - (b./r).^2 - (g^-2)*(2/15*(s/R)^9 *(1./(r - 1).^9 - 1./(r + 1).^9 - 9./(8* r).*(1./(r - 1).^8 - 1./(r + 1).^8)) - (s/R)^3 *(1./(r -1).^3 - 1./(r + 1).^3 - 3./(2* r).* (1./(r - 1).^2 - 1./(r + 1).^2)))));
xx = real ( pi - 2*b * integral(fun,max(vec),Inf));

end
and it evaluates very quickly de double integral, but when I started to evaluate the triple integral, that is an integral of a double integral qQd, it doesnt work, take a lot of time

% integral colisao difusa (1,1) 

function exem8  
tic

T=100;

   f3 =@(g)  g.^5.*qQd2(g)./(exp(g.^2/T));
  od11= (1/T^3)*integral(f3,0,1)
  toc
      end

      function Res2 = qQd2(g)
       for i=1:length(g)
          Res2(i)=qQd(g(i));
      end
      end

      function Res = qQd(g)

      f=@(b) anggbsR(b,g);
  %x0=[1,2];
   x=fzero(f,1.1);
  %T=100;
  %g=10;

     T=100; 
      f1=@(b) 2*(1 + (sqrt(pi*T)/(2*g)).*sin(ang2(b,g)/2)).*b ;
  f2=@(b)  2*(1 - cos(ang2(b,g))).*b ;
  Res = integral (f2,x, 10,'RelTol',0,'AbsTol',1e-7) + integral(f1,0,x);
      end

      function xx2 = ang2(b,g)
      for i=1:length(b)
          xx2(i)=anggbsR(b(i),g);
       end
      end

function xx=anggbsR(b,g)

T=100;
s= 0.6; 
R= 1; 
syms r
func =  @(r)  1 - (b/r)^2 - (g^-2)*(2/15*[(s/R)]^9 *(1/(r - 1)^9 - 1/(r + 1)^9 - 9/(8* r)*(1/(r - 1)^8 - 1/(r + 1)^8)) - [(s/R)]^3 *(1/(r -1)^3 - 1/(r + 1)^3 - 3/(2* r)* (1/(r - 1)^2 - 1/(r + 1)^2)));
minrootgbsR = fzero(func,[1.01,100]);
k=1;
for i=1:length(minrootgbsR)
    if(isreal(minrootgbsR(i))==1)
        vec(k)=minrootgbsR(i);
        k=k+1;
    end

end

fun =@(r) 1./(r.^2 .* sqrt(1 - (b./r).^2 - (g^-2)*(2/15*(s/R)^9 *(1./(r - 1).^9 - 1./(r + 1).^9 - 9./(8* r).*(1./(r - 1).^8 - 1./(r + 1).^8)) - (s/R)^3 *(1./(r -1).^3 - 1./(r + 1).^3 - 3./(2* r).* (1./(r - 1).^2 - 1./(r + 1).^2))))); 
xx = real ( pi - 2*b * integral(fun,max(vec),Inf));

end
and give this errors :

Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   1.0e-05. The integral may not exist, or it may be difficult to approximate numerically to the
requested accuracy. 
> In funfun\private\integralCalc>iterateScalarValued at 372
  In funfun\private\integralCalc>vadapt at 132
  In funfun\private\integralCalc at 83
  In integral at 88
  In exem8>anggbsR at 66
  In exem8>ang2 at 38
  In exem8>@(b)2*(1-cos(ang2(b,g))).*b at 31
  In funfun\private\integralCalc>iterateScalarValued at 323
  In funfun\private\integralCalc>vadapt at 132
  In funfun\private\integralCalc at 75
  In integral at 88
  In exem8>qQd at 32
  In exem8>qQd2 at 15
  In exem8>@(g)g.^5.*qQd2(g)./(exp(g.^2/T)) at 8
  In funfun\private\integralCalc>iterateScalarValued at 314
  In funfun\private\integralCalc>vadapt at 132
  In funfun\private\integralCalc at 75
  In integral at 88
  In exem8 at 9

[Walter Roberson]

You should probably be using integral2() and integral3(). You can pass function handles for the bounds for the second and third variables to handle the root finding you are doing.