integration for spectral radiance formula I have the following formula:
$$ L_{\lambda} = \frac{2c^{2}h}{\lambda^{5}\left( \exp\left( \frac{hc}{\lambda kT} \right)-1 \right)}$$
which is used to calculate spectral radiance (see here).
I have a matlab function for calculating the spectral radiance across all wavelengths:
function b6000 = myfun(lambda) 

h = 6.626e-34; % Planck's Constant
c = 3e8; % speed of light
T = 6000;
k = 1.38066e-23; % Boltzmann constant in J/K

% spectral radiance
p = 2*h*c*c./(lambda.^5);
b6000 = p./(exp((h*c)./(lambda*k*T))-1);
b6000 = (1e-9).*b6000;

% multiply by the square of the ratio of the solar radius of earth's
% orbital radius
b6000 = b6000.*(2.177e-5);

% apply Lambert's cosine law
b6000 = b6000.*pi;

where the entire area under the spetral radiance curve can be computed by:
integral(@myfun,0,inf)

which gives an answer of 
1.5977e-06

which is in the wrong units to make any sense, but is the correct answer. 
I am trying to calculate this by hand, but am struggling to know where I would start.
$$ L_{\lambda} = \int_{0}^{\inf} \frac{2c^{2}h}{\lambda^{5}\left( \exp\left( \frac{hc}{\lambda kT} \right)-1 \right)} d\lambda$$
I do not know of any methods than can be used to integrate this. I am not a mathematician, however, so there most likely is a sensible way of doing this. Can someone provide any pointers? 
 A: You can basically use Gaussian quadrature. Or the command "quad" on Matlab.
A: Basically the integral is (ignoring the constant terms) 
$$I=\int_0^{\infty} \frac{d\lambda}{\lambda^5 \left(\exp\left(\frac{a}{\lambda}\right)-1\right)}$$
where $a=\dfrac{hc}{kT}$. Substitute $t=a/\lambda \Rightarrow d\lambda=-a\,dt/t^2$ to obtain the integral:
$$I=\frac{1}{a^4}\int_0^{\infty} \frac{t^3}{e^t-1}\,dt=\frac{1}{a^4}\int_0^{\infty} \frac{t^3e^{-t}}{1-e^{-t}}\,dt$$
Next use the series expansion $\displaystyle \frac{1}{1-e^{-t}}=\sum_{k=0}^{\infty} e^{-kt}$ to rewrite the integral as:
$$I=\frac{1}{a^4}\sum_{k=0}^{\infty} \int_0^{\infty} t^3e^{-(k+1)t}\,dt$$
Use the substitution $(k+1)t=x$,
$$\Rightarrow I=\frac{1}{a^4}\sum_{k=0}^{\infty} \frac{1}{(k+1)^4}\int_0^{\infty}x^3e^{-x}=\frac{1}{a^4}\sum_{k=0}^{\infty} \frac{1}{(k+1)^4}\Gamma(4)=\frac{3!}{a^4}\sum_{k=0}^{\infty} \frac{1}{(k+1)^4}$$
$$\Rightarrow I=\frac{6}{a^4}\sum_{k=1}^{\infty} \frac{1}{k^4}\frac{6}{a^4}\frac{\pi^4}{90}=\frac{\pi^4}{15a^4}$$
Hence, the final answer is:
$$2hc^2I=2hc^2\frac{\pi^4}{15}\frac{k^4T^4}{h^4c^4}=\boxed{\dfrac{2\pi^4k^4}{15h^3c^2}T^4}$$
