Planck's law question I was reading Wikipedia article about Planck's Law and I wanted to make the same graph as here.
I took this equation
$$ B(\nu,T)=\frac{2h\nu^3}{c^2}\frac{1}{e^{\frac{h\nu}{kT}}-1}$$
but I got different graphs. Maximum is shifted. For example, in Wikipedia maximum for 5000K is around 0.5-0.6$\mu m$, but I got maximum for 5000K at 1$\mu m$. My graph is below:

I used this Wolfram Mathematic program below to build this graph.
kV = 1.380649*10^(-23);
hV = 6.62607015*10^(-34);
cV = 299792458;
f[l_] = cV/l/10^(-6);(*Convert Wave Length to Frequency*)
Bf[f_, T_] = 2*hV*f^3/cV^2*1/(Exp[(hV*f)/(kV*T)] - 1);

Plot[{Bf[f[l], 5000], Bf[f[l], 4000], Bf[f[l], 3000]},
 {l, 0.3, 2},
 PlotStyle -> {Blue, Green, Red},
 PlotLegends -> {"5000K", "4000K", "3000K"},
 ImageSize -> 800]

My question is what am I doing wrong and why I am getting a different maximum?
UPDATE 1.
I removed $\mu m$ conversion to have everything in the same SI units.
Graph and program below.

kV = 1.380649*10^(-23);
hV = 6.62607015*10^(-34);
cV = 299792458;
f[l_] = cV/l;
Bf[f_, T_] = 2*hV*f^3/cV^2*1/(Exp[(hV*f)/(kV*T)] - 1)

Plot[{Bf[f[l], 5000], Bf[f[l], 4000], Bf[f[l], 3000]},
 {l, 0.3*10^(-6), 2*10^(-6)},
 PlotStyle -> {Blue, Green, Red},
 PlotLegends -> {"5000K", "4000K", "3000K"},
 ImageSize -> 800]

 A: Ah you got caught by something subtle! Planck's law gives a probability density function in terms of the frequency $B(\nu,T)$. To convert to the probability density function in terms of $\lambda$ (let's call it $\tilde{B}(\lambda,T)$), you can't just set $\nu=c/\lambda$. You need to include a Jacobian factor.
In more detail:
What remains invariant under a change in variables is the probability $dp=B(\nu,T)d\nu=\tilde{B}(\lambda,T)d\lambda$. As a result...
\begin{equation}
\tilde{B}(\lambda,T) = \left|\frac{d \nu}{d\lambda}\right| B(\nu,T) = \frac{c}{\lambda^2} B\left(\frac{c}{\lambda},T\right)
\end{equation}
I bet if you include the $c/\lambda^2$ factor in your plot, you will get the same result as wikipedia. Let me know how it goes!
A: You are comparing the spectral radiance in frequency $B(f,T)$ with the spectral radiance in wavelength $B(\lambda,T)$. Because these are differential distributions, you cannot simply replace $f=c/\lambda$ to go from one to another. You must take into account the Jacobian of transformation.
$$d\Omega=B(f,T)df = B(\lambda,T)d\lambda$$
$$\begin{align}
\implies B(\lambda,T)&=\overset{\color{red}{[1]}}{-}B(f(\lambda),T)\frac{df}{d\lambda}\\
&=+\frac{c}{\lambda^2}B(f(\lambda),T)\\
&=+\frac{f^2}{c}B(f(\lambda),T)
\end{align}$$
If you change that definition it should work fine.
Edit: The $\color{red}{[1]}$ marker above is in reference to @rob's comment, which is a minus sign I forgot to write. One must take into account limits of integration in the above relation between $B(\lambda,T)$ and $B(f,T)$. If $f_1<f_2$, then $\lambda_1>\lambda_2$, and therefore:
$$\begin{align}
\int_{\lambda_2}^{\lambda_1} B(\lambda,T)d\lambda &= \int_{f_1}^{f_2} B(f,T)df\\
&=\int_{\lambda(f_1)=\lambda_1}^{\lambda(f_2)=\lambda_2} B(f(\lambda),T)
\frac{df}{d\lambda}\,d\lambda \\
&=\int_{\lambda_2}^{\lambda_1} -B(f(\lambda),T)
\frac{df}{d\lambda}\,d\lambda \\
&\\
& \implies \boxed{B(\lambda,T)=-B(f(\lambda),T)\frac{df}{d\lambda}}
\end{align} $$
