Solving coupled 2nd order ODE, numerical in Python

I would like to solve the following DGL system numerically in python:

The procedure should be the same as always. I have 2 coupled differential equations of the 2nd order and I use the substitution g' = v and f' = u to create four coupled differential equations of the 1st order.

Here ist my code:

from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt

# Constants
e = 1
mu = 1
lambd = 1

# Initial second-order system

# g'' + 2/r g'-2/r^2 g - 3/r eg^2- e^2 g^3 - e/(2r) f^2 - e^2/4 f^2g = 0
# f'' + 2/r f' - 2/r^2 f - 2e/r fg - e^2/2 fg^2 + mu^2 f - lambda f^3 = 0

# Substitution
# g' = v
# v' = g''
# f' = u
# u' = f''

# Equivalent set of first-order systems
def dSdr(r, S):
    g, v, f, u = S
    dgdr = v
    dvdr = -2 / r * v + 2 / r ** 2 * g + 3 / r * e * g ** 2 + \
        e ** 2 * g ** 3 + e / (2 * r) * f**2 + e**2 / 4 * f ** 2 * g
    dfdr = u
    dudr = -2 / r * u + 2 / r ** 2 * f + 2 * e / r * f * g + e ** 2 / 2 * f * g**2 - mu ** 2 * f + lambd * f ** 3

    return [dgdr, dvdr, dfdr, dudr]

# Initial conditions
g0 = 0.0001
v0 = 0.0001
f0 = 1
u0 = 0.0001
S0 = [g0, v0, f0, u0]

r_start = 0.1
r_end = 100
r = np.linspace(r_start, r_end, 1000)

# Solve the differential equations

sol = solve_ivp(dSdr, t_span=(min(r), max(r)), y0=S0, t_eval=r, method='RK45')




# Check if integration was successful
if sol.success:
    print("Integration successful")
else:
    print("Integration failed:", sol.message)

# Debugging information
print("Number of function evaluations:", sol.nfev)
print("Number of accepted steps:", sol.t_events)
print("Number of steps that led to errors:", sol.t)

# Plot results
plt.figure(figsize=(10, 6))
plt.plot(sol.t, sol.y[0], label='g(r)')
plt.plot(sol.t, sol.y[2], label='f(r)')
plt.xlabel('r')
plt.ylabel('Function values')
plt.legend()
plt.title('Solutions of the differential equations')
plt.show()

My boundary con. should be f(0) = g(0) = 0 and f(infinity) = mu/sqrt(lambd) g(infinity) = 0 so that the system makes physical sense. But how can I incorporate this condition or how do I know the initial conditions for v,u? The system should look like this (From the original paper):

but it doesn't. Does anyone know what to do?

Source:

enter link description here

Although you could try a "shooting" method using initial-value solvers, I find it best to go for a boundary-value problem directly. I do so by a method akin to the finite-volume method in computational fluid mechanics (CFD).

Division by r creates singularities, so start by multiplying your equations by r². Collecting the derivative terms you will then get equations which look like spherically-symmetric diffusion equations:

Here, S_f and S_g are the great long source terms (which you will probably have to check in my code - there's plenty of scope for minor errors). Discretise each (equivalent to integrating over a cell [r_i-1/2,r_i+1/2] by the finite-volume method). e.g. for f:

This, and the corresponding equation for g can be written in the form

There are small differences at the boundaries and it is the non-zero boundary condition on f that produces a non-trivial solution.

The discretised equations can then be solved iteratively either by the tri-diagonal matrix algorithm (TDMA) or by iterative schemes like Gauss-Seidel or Jacobi. Although experience with CFD suggests that the TDMA would probably be best, for simplicity I’ve opted for a (slightly under-relaxed) Jacobi method here.

Note that, in common with CFD practice, I have, in my edited code, split the source terms as, e.g. sf + spf.f, with spf being negative so that it can be combined with the diagonal coefficient ap on the LHS.

EDIT: code amended so that cell-centre values are 1...N with boundary indices 0 and N+1, source terms split into explicit and implicit parts, code reads more like Python than Fortran!

import numpy as np
import matplotlib.pyplot as plt

N = 100                                # number of cells (1-N); added boundaries 0 and N+1
rmin, rmax = 0.0, 20.0                 # minimum and maximum r (latter an approximation to "infinity"
dr = rmax / N                          # cell width
gL, gR, fL, fR = 0, 0, 0, 1            # boundary conditions

e, nu, lamda = 1, 1, 1

# Cell / node arrangement:
#
#       0                                          N+1
#       |  1     2                        N-1    N  |
#       |-----|-----| .  .  .  .  .   . |-----|-----|
#
r = np.linspace( rmin - 0.5 * dr, rmax + 0.5 * dr, N + 2 )      # cell centres (except for boundaries)
r[0] = rmin;   r[N+1] = rmax                                    # boundary nodes on faces of cells

# Connection coefficients
awL = 2 * rmin ** 2 / dr ** 2
aeR = 2 * rmax ** 2 / dr ** 2
aw = np.zeros( N + 2 )
ae = np.zeros( N + 2 )
for i in range( 1, N + 1 ):
    if i == 1:
        aw[i] = 0
    else:
        aw[i] = ( 0.5 * ( r[i-1] + r[i] ) ) ** 2 / dr ** 2
    if i == N:
        ae[i] = 0
    else:
        ae[i] = ( 0.5 * ( r[i+1] + r[i] ) ) ** 2 / dr ** 2
ap = aw + ae

# Initial values
g = np.zeros( N + 2 )
f = np.zeros( N + 2 )
f = f + 1
g = g + 1

# Boundary conditions
f[0] = fL;   f[N+1] = fR;   g[0] = gL;   g[N+1] = gR


alpha = 0.9                  # Under-relaxation factor
niter = 0
for _ in range( 10000 ):
    niter += 1

    # Source term (e.g., for f this would be decomposed as sf + spf.f, where spf is negative)
    spf = -2  -0.5 * e**2 * r**2 * g**2  - lamda * r**2 * f**2
    sf = -2 * e * r * f * g   +   nu**2 * r**2 * f

    spg = -2  - e**2 * r**2 * g**2  - 0.25 * e**2 * r**2 * f**2
    sg  = -e * r * ( 3 * g**2 + 0.5 * f**2 )

    # Dirichlet boundary conditions applied via source term
    sf[1] += awL * fL;   spf[1] -= awL
    sg[1] += awL * gL;   spg[1] -= awL
    sf[N] += aeR * fR;   spf[N] -= aeR
    sg[N] += aeR * gR;   spg[N] -= aeR

    # Update g and f (under-relaxed Jacobi method)
    ftarget = f.copy()
    gtarget = g.copy()
    for i in range( 2, N ):                                # cells 2 - (N-1)
        ftarget[i] = ( sf[i] + aw[i] * f[i-1] + ae[i] * f[i+1] ) / ( ap[i] - spf[i] )
        gtarget[i] = ( sg[i] + aw[i] * g[i-1] + ae[i] * g[i+1] ) / ( ap[i] - spg[i] )
    i = 1                                                  # leftmost cell
    ftarget[i] = ( sf[i] + ae[i] * f[i+1] ) / ( ap[i] - spf[i] )
    gtarget[i] = ( sg[i] + ae[i] * g[i+1] ) / ( ap[i] - spg[i] )
    i = N                                                  # rightmost cell
    ftarget[i] = ( sf[i] + aw[i] * f[i-1] ) / ( ap[i] - spf[i] )
    gtarget[i] = ( sg[i] + aw[i] * g[i-1] ) / ( ap[i] - spg[i] )

    # Under-relax the update
    f = alpha * ftarget + ( 1 - alpha ) * f
    g = alpha * gtarget + ( 1 - alpha ) * g

    change = ( np.linalg.norm( f - ftarget ) + np.linalg.norm( g - gtarget ) ) / N
    # print( niter, change )
    if change < 1.0e-10: break

print( niter, " iterations " )

plt.plot( r, f, label='f' )
plt.plot( r, g, label='g' )
plt.grid()
plt.legend()
plt.show()