Solution program [deuteron.f90 (Fortran)] [Julia code (Jupyter notebook) by Gabe Schumm)]
The results shown below were obtained using rmax=100 fm, and 20000 integration steps, but even rmax=50 and 2000 steps give acceptable results.