Help understanding why quadgk/quadl does not work

Im trying to integrate a specific function, but it is giving me a hard time getting it to work. I want to integrate a product of functions. Im using the build in Regularized Incomplete Beta Function betainc and an own function to calculate the PDF/CDF of a Generalized Hyperbolic Distribution. For completeness of my question, this is the code for the PDF:

function y=ghyppdf(t,lam,alpha,beta,delta,mu)
%Generalized Hyperbolic probability density function (pdf).
%   Y=GHYPPDF(X,LAM,ALPHA,BETA,DELTA,MU) returns the pdf of the generalized
%   hyperbolic distribution with parameter LAM,
%   shape parameter ALPHA, skewness parameter BETA, 
%   scale parameter DELTA and location parameter MU, evaluated at the 
%   values in X.

if (delta>0)&(alpha>abs(beta))
    kappa = (alpha^2-beta^2)^(lam/2)/(sqrt(2*pi)*alpha^(lam-1/2)*delta^lam*besselk(lam,delta*sqrt(alpha^2-beta^2)));
    y = kappa*(delta^2+(t-mu).^2).^((lam-1/2)/2).*besselk(lam-1/2,alpha*sqrt(delta^2+(t-mu).^2)).*exp(beta*(t-mu));
else
    y = Inf*ones(size(t));
end;
y(isnan(y))=0;
end

And the CDF:

function y=ghypcdf(x,lam,alpha,beta,delta,mu,starti)
%HYPCDF Hyperbolic cumulative distribution function (cdf).
%   Y=HYPCDF(X,ALPHA,BETA,DELTA,MU,STARTI) returns the cdf of the hyperbolic  
%   distribution with shape parameter ALPHA, skewness parameter BETA, 
%   scale parameter DELTA and location parameter MU, evaluated at the 
%   values in X. STARTI is an optional parameter specifying the starting
%   point for the integration scheme.

% Convert to a column vector
x = x(:);

% Find the starting point for the integration scheme
feps = 1e-10;
if nargin < 7
    starti = mu+beta*delta/sqrt(alpha^2-beta^2)*besselk(lam+1,delta*sqrt(alpha^2-beta^2))/besselk(lam,delta*sqrt(alpha^2-beta^2));
    starti = min(starti,min(x))-1;
    while ghyppdf(starti,lam,alpha,beta,delta,mu)>feps
        starti = starti-1;
    end;
end;

n = length(x);
y = zeros(n,1);
[x,ind] = sort(x);
ind = sortrows([ind,(1:n)'],1);
ind = ind(:,2);
x = [starti;x];
warning off MATLAB:quadl:MinStepSize
% Integrate the hyperbolic pdf
for i=1:n
    y(i) = quadl('ghyppdf',x(i),x(i+1),[],[],lam,alpha,beta,delta,mu);
end;
warning on MATLAB:quadl:MinStepSize
y = cumsum(y);
y = y(ind);
y(y<0) = 0;%security
y(y>1) = 1;
y(isnan(y)) = 0;%security
end

Now im trying to calculate the following integral in three ways:

1) EV = integral( @(a)betainc(1-ghypcdf(a,lam,alpha,beta,delta,mu),T-Tp,Tp).*a.*ghyppdf(a,lam,alpha,beta,delta,mu),-Inf,Inf,'ArrayValued',true)

2) EV = quadgk( @(a)betainc(1-ghypcdf(a,lam,alpha,beta,delta,mu),T-Tp,Tp).*a.*ghyppdf(a,lam,alpha,beta,delta,mu),-Inf,Inf)

3) EV = quadl( @(a)betainc(1-ghypcdf(a,lam,alpha,beta,delta,mu),T-Tp,Tp).*a.*ghyppdf(a,lam,alpha,beta,delta,mu),-Inf,Inf)

For some reason only the first method works while the other two give me errors saying that the matrix dimensions do not agree. Altough method 1 is working, i want to understand what is going wrong with the others. Based on this i have two questions:

1) Why do i need to include the 'ArrayValued',true part in option 1? It looks like matlab is treating something as an vector/matrix, but it is clearly a scalar function?

2) Why are method 2 and 3 not working, while i explicitly use .* to make a matrix product?

Thanks in advance!