The system (1)-(5) was integrated, using a semi-implicit Euler
method. A uniform rectangular grid with radial and angular steps 
in a rectangle 
was used. The derivatives with respect to r
and 
were approximated by central differences, except for the
radial derivative 
at the boundaries, where the
one-sided second order differences were used
To provide convergence at small r, the simplest implicit finite difference approximation of the angular part of the Laplacian was used. Both the radial part of the Laplacian and the non-linear terms were integrated explicitly. To resolve the resultant linear system a cyclic elimination method was implemented [24]. The Laplacian at r=0 was computed by formula (8) resulting from the mean value theorem for a disk [11], [22], [23],
More precisely, the non-local boundary condition at r=0 has been set
Here 
are the values of u at the circumference with radius
, M is the number of angular segments, and 
is the value
of u at r=0.
