How to solve NaN values of integrand defined over a finite interval?

Hi, could someone know why quadv () returns NaN (see %r1), although my defined function is well behaved on [-theta; pi/2]?

see below the code:

f=@(k) k .* exp(-k);
x=1
al=1;
bet=0.5;

if al==1
    c=exp(-(pi * x) /2 * bet);
    c2=1 /(2*abs(bet));
    theta0=pi/2
    zet=0
    v1=@(theta) 2/pi .* ((pi/2 + bet.*theta) ./cos(theta)) .* exp(1/bet .* (pi/2 + bet.*theta)...
    .* tan(theta));
    g=@(theta) c .* v1(theta);
    %p=c2 * quadv(@(theta) f(g(theta)),-theta0,pi/2);
    y=  @(theta) c2 .* f(g(theta))    ;  % values of integrand
    %r1= quadv(@(theta) y(theta),-theta0,pi/2)

      % check behavior of integrand

      n= @(theta) f(g(theta));
      a=[-theta0:0.01:pi/2]';
      N=n(a);
      P=[a;N];
      fprintf('%6s %12s\n','x','f(x)');
      fprintf('%6.2f %12.8f\n',P);
      plot(a,P)

else
    alm=al-1;
    d=1/alm;
    da=al .*d;
    theta0=(1/al) .*atan(bet .*tan(pi .*al/2));
    u=al .*theta0;
    zet=-bet .*tan(pi*al/2)
    v=@(theta) (cos(u) .^d) .*((cos(theta) ./sin(al .*(theta0+theta))) .^da) .*...
        (cos(u+alm .*theta) ./cos(theta));
    c=(x-zet) .^da;
    c2=al ./(pi .*abs(al-1) .* (x-zet));
    g=@(theta) c .* v(theta);
    y=  @(theta) f(g(theta))    ;  % values of integrand
    %y=  @(theta) c2 .* f(g(theta))    ;
    r2= c2 .* quadv(@(theta) y(theta),-theta0,pi/2)
end

% plot values
%a=[-theta0:0.01:pi/2]';
%Y=y(a);
%U=[a;Y];
%fprintf('%6s %12s\n','x','f(x)');
%fprintf('%6.2f %12.8f\n',U);
%plot(a,Y)

Also for curiosity I am wondering why from the screen x goes beyond the maximum limit pi/2. actually matlab outputs lots of zeros there.

Thanks