  The closed-form solution of the 1.5 post-Newtonian (PN) accurate binary black
hole (BBH) Hamiltonian system has proven to be difficult to obtain for a long
time since its introduction in 1966. Closed-form solutions of the PN BBH
systems with arbitrary parameters (masses, spins, eccentricity) are required
for modeling the gravitational waves (GWs) emitted by them. Accurate models of
GWs are crucial for their detection by LIGO/Virgo and LISA. Only recently, two
solution methods for solving the BBH dynamics were proposed in arXiv:1908.02927
(without using action-angle variables), and arXiv:2012.06586, arXiv:2110.15351
(action-angle based). This paper combines the ideas laid out in the above
articles, fills the missing gaps and provides the two solutions which are fully
1.5PN accurate. We also present a public Mathematica package BBHpnToolkit which
implements these two solutions and compares them with a fully numerical
treatment. The level of agreement between these solutions provides a numerical
verification for all the five actions constructed in arXiv:2012.06586, and
arXiv:2110.15351. This paper hence serves as a stepping stone for pushing the
action-angle-based solution to 2PN order via canonical perturbation theory.
