vpa() is for use with symbolic expressions only.
If you have the symbolic engine,
besselFunction = @(u)besseli(qb,2*sqrt(lambda*(theta + mu)).*sym(u));
exponentFuncion = @(u)exp(-sym(u).*(lambda + theta + mu));
Then besselFunction(u)*exponentFunction(u) will return a symbolic value. You can double() it if you want the double-precision representation.
Note: going symbolic will only postpone the problem: the symbolic engine is limited to about a billion decimal places.
Best Answer