  Systems of non-zero isospin chemical potential are studied from a canonical
approach by computing correlation functions with the quantum numbers of $N
\pi^+$'s ($C_{N \pi}$). In order to reduce the number of contractions required
in calculating $C_{N \pi}$ for a large $N$ in the Wick's theorem, we
constructed a few new algorithms. With these new algorithms, systems with
isospin charge up to 72 are investigated on three anisotropic gauge ensembles
with a pion mass of $390 \rm{MeV}$, and with lattice spatial extents $L \sim
{2.0, 2.5, 3.0} \rm{fm}$. The largest isospin density of $\rho_I \approx 9
\rm{fm}^{-3}$ is achieved in the smallest volume, and the QCD phase diagram is
investigated at a fixed low temperature at varying isospin chemical potentials,
$m_{\pi} \le \mu_I \le 4.5 m_{\pi}$. By investigating the behaviour of the
extracted energy density of the system at different isospin chemical
potentials, we numerically identified the conjectured transition to a
Bose-Einstein condensation state at $\mu_I \ge m_{\pi}$.
