% 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}
