Harmonic oscillator energy vs inverse temperature beta = 1/(kB T)
\documentclass{standalone}
\usepackage{pgfplots}
\pgfplotsset{compat=newest}
\begin{document}
\begin{tikzpicture}
\begin{axis}[
xlabel = $\beta$,
ylabel = $\langle E\rangle/\hbar \omega$,
ymin = 0,ymax = 2.3,
domain = 0:11,
smooth,thick,
axis lines = center,
every tick/.style = thick]
\def\h{1}\def\w{1}
\def\eev{1/2*\h*\w*(1 + 4/(e^(x*\h*\w) - 1))}
\addplot[color=blue]{\eev};
\end{axis}
\end{tikzpicture}
\end{document}