Wu-Yang monopole: numerical calculation

I have been reading the paper by Wu and Yang [1] in which they find the famous Wu-Yang monopole. In the paper there are solutions for three types of monopoles: one has an analytical form, which is the one most often quoted, but there are also two other monopoles with numerical solution only. In this post I use Python/numpy to perform numerical analysis on the latter solution. I use the same notation as in [1].
Wu and Yang obtain the following system of ordinary differential equations $$\begin{align} \frac{d\Phi}{d \xi} &= \psi\label{eq:20170625a}\\ \frac{d\psi}{d \xi} &= \psi + \Phi(\Phi^2-1)\label{eq:20170626a} \end{align}$$ Here $\xi$ is given by $r = e^{\xi}$, with $r$ the distance to the origin. The right-hand side of $\eqref{eq:20170625a}$-$\eqref{eq:20170626a}$ defines the vector field ($d\Phi/d\xi, d\psi/d\xi)$ in the $(\Phi, \psi)$ plane. Its integral curves can be seen in the next figure

The integral curves of the vector field defined by $\eqref{eq:20170625a}$-$\eqref{eq:20170626a}$.
The stationary points are marked in red.

I calculate the integral curve from the point $(\Phi,\psi) = (0,0)$ to $(1,0)$ using the numpy function solve_bvp [2].

The integral curves of the vector field defined by $\eqref{eq:20170625a}$-$\eqref{eq:20170626a}$.
The integral curve from the stationary point $(0,0)$ to $(1,0)$ is added in red.

$\Phi(\xi)$ can be seen in the next graph. One sees that $\Phi(\xi) \to 0$ for $\xi \to -\infty$ and $\Phi(\xi) \to 1$ for $\xi \to +\infty$

In the rest of this post I reproduce part of Table 1 in [1].
Notice that if $(\Phi(\xi),\psi(\xi) )$ is a solution of $\eqref{eq:20170625a}$-$\eqref{eq:20170626a}$, then $(\Phi(\xi- \xi_0),\psi(\xi-\xi_0) )$ with $\xi_0$ a constant is also a solution. Translating to the $r$ coordinate, this means that if $\Phi(r)$ is a solution, the rescaled function $\Phi(r/c)$ is also a solution. Wu and Yang provide a table with numerical results on the function $\Phi(r)$ with asymptotic behaviour $\Phi(r) = 1 - 1/r +O(1/r^2)$ for $r \to \infty$. The solution plotted above has $\Phi(r) = 1 - c/r+O(1/r^2)$ for $r \to \infty$ with $c = 0.6233$ [3]. If I rescale my solution with $c$, I get the following table.

$\xi$	$r/c$	$\Phi(r/c)$
5.066	9.880e+01	9.898e-01
2.866	1.095e+01	9.136e-01
1.666	3.297e+00	7.508e-01
0.566	1.098e+00	4.583e-01
-2.317	6.141e-02	-9.229e-02
-5.954	1.617e-03	1.498e-02
-9.583	4.296e-05	-2.442e-03

This agrees well with the results in Table 1 of Wu and Yang. I have not compared smaller values of $r$ because I have approximated the infinite interval $-\infty < \xi < + \infty$ by $-12\le \xi \le 8$.

References and comments

[1] Wu and Yang, Some Solutions of the Classical Isotopic Gauge Field Equations, 1969

[2] The Python code that I used can be found at this link.

[3] I have not estimated $c$ from the asymptotic behaviour of the solution that I found, but I have taken $c$ so that the difference between my solution and the one in Wu Yang is as small as possible.