The values are so small relative to your constants that they might as well be zero.
With your beta being 1/(k*t) you are ending up with some large exponents in the exp() -- large enough that you are going to overflow numerically.
syms x
Q = @(v) sym(v);
k = Q(1380649)*Q(10)^(-29);
t = Q(5800);
beta = (k*t)^(-1);
fun = x.^3 ./ (exp(beta .* x) - 1)
int(fun, 0, Inf)
vpa(ans)
