Quarkonium is an ideal probe to explore the properties of quantum chromodynamics (QCD). Unlike Large Hadron Collider (LHC) measurements, quarkonium production at the Relativistic Heavy Ion Collider (RHIC) has different production mechanisms, can access different kinematic phase space and may experience different medium densities/temperatures. The PHENIX experiment has collected a large $J/\psi \to \mu^{+}\mu^{-}$ data set within its unique pseudorapidity region of $1.2<|\eta|<2.2$ in $p+p$, $p+A$ and $A+A$ collisions at $\sqrt{s}$ = 200 GeV from 2014 to 2016. Latest results of normalized charged particle multiplicity dependent normalized forward $J/\psi$ yields and normalized forward $\psi(2S)$ to $J/\psi$ ratio in 200 GeV $p+p$ collisions as well as forward $J/\psi$ azimuthal anisotropy in 200 GeV Au+Au collisions will be shown. Comparison with other RHIC and LHC measurements and latest theoretical calculations will be discussed. These PHENIX results provide their unique contributions in improving the understanding of the multi-parton interaction contribution to charmonium production and the recombination/coalescence contribution to charmonium formation within Quark Gluon Plasma (QGP) at RHIC energies.