How to evaluating solutions to Coupled Nonlinear Elliptic ODEs with Relaxation Method?

I would like to determine the following system of coupled, nonlinear, elliptic ODEs of second order with boundary conditions f(0) = h(0) = 0, f(1) = h(1) = 1.

I used a relaxation method to solve the system (at least that's what I understand by it). However, my code gives me a rather caotic result, which cannot be correct. I am not used to dealing with such difficult differential equations, is another method perhaps more suitable? Perhaps someone sees a mistake that I have made.

import numpy as np
import matplotlib.pyplot as plt

# Define parameters
N = 100  # Number of points
lambda_val = 1  # Example value for lambda
g = 1  # Example value for g
tolerance = 1e-6  # Convergence tolerance
max_iterations = 10000  # Maximum number of iterations

# Discretize the domain
x = np.linspace(0, 1, N + 1)
dx = x[1] - x[0]

# Initial guesses
f = x.copy()
h = x.copy()

# Boundary conditions
f[0], h[0] = 0, 0
f[1], h[1] = 1, 1

# Relaxation method
for iteration in range(max_iterations):
    f_old = f.copy()
    h_old = h.copy()

    # Update f
    for i in range(1, N):
        # Discretized derivatives
        f_prime = (f[i + 1] - f[i - 1]) / (2 * dx)
        f_double_prime = (f[i + 1] - 2 * f[i] + f[i - 1]) / (dx ** 2)

        # Avoid division by zero
        h_term = x[i] ** 2 * h[i] ** 2 / (x[i] - 1) ** 2 if x[i] != 1 else 0
        rhs_f = (1 - f[i]) * (8 * (1 - 2 * f[i]) * f[i] - h_term) / 4
        lhs_f = 2 * (x[i] - 1) * x[i] ** 2 * f_prime + (x[i] - 1) ** 2 * x[i] ** 2 * f_double_prime

        # Ensure the update step is reasonable
        update_f = (rhs_f - lhs_f) * dx ** 2
        if np.abs(update_f) < 1e6:  # Avoid too large updates
            f[i] = f_old[i] + update_f

    # Update h
    for i in range(1, N):
        # Discretized derivatives
        h_prime = (h[i + 1] - h[i - 1]) / (2 * dx)
        h_double_prime = (h[i + 1] - 2 * h[i] + h[i - 1]) / (dx ** 2)

        # Avoid division by zero
        f_term = 2 * (f[i] - 1) ** 2
        lambda_term = x[i] ** 2 * lambda_val * (h[i] ** 2 - 1) / (g ** 2 * (x[i] - 1) ** 2) if x[i] != 1 else 0
        rhs_h = h[i] * (f_term + lambda_term)
        lhs_h = 2 * (x[i] - 1) ** 2 * x[i] * h_prime + (x[i] - 1) ** 2 * x[i] ** 2 * h_double_prime

        # Ensure the update step is reasonable
        update_h = (rhs_h - lhs_h) * dx ** 2
        if np.abs(update_h) < 1e6:  # Avoid too large updates
            h[i] = h_old[i] + update_h

    # Check for convergence
    if np.max(np.abs(f - f_old)) < tolerance and np.max(np.abs(h - h_old)) < tolerance:
        print(f"Converged after {iteration} iterations.")
        break
else:
    print("Did not converge within the maximum number of iterations.")

# Print the final solution
print("f =", f)
print("h =", h)

# Plotting the results
plt.figure(figsize=(12, 6))

# Plot f(x)
plt.subplot(1, 2, 1)
plt.plot(x, f, label='f(x)')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Solution for f(x)')
plt.legend()
plt.grid(True)

# Plot h(x)
plt.subplot(1, 2, 2)
plt.plot(x, h, label='h(x)', color='orange')
plt.xlabel('x')
plt.ylabel('h(x)')
plt.title('Solution for h(x)')
plt.legend()
plt.grid(True)

plt.tight_layout()
plt.show()

As I put in a comment, I don’t think you should do a transformation for x just to bring infinity to a finite point: I think you should stick with the original equations in r (equations 11.44 and 11.45 in your linked document) and just put the right boundary at a reasonably large number. I have solved these original equations by the finite-volume method in a slightly improved version of that for one of your previous threads, Solving coupled 2nd order ODE, numerical in Python

Slightly rearranged – for reasons given below – the original equations are

These CONSERVATIVE equations can be discretised by integrating over a cell. Note that I divide the source term on the RHS into an explicit and an implicit part; the coefficient for the latter must be negative.

Each of these can then be written as a TRI-DIAGONAL SYSTEM:

where

Making the implicit source term negative ensures that this is diagonally dominant. The boundary conditions are usually handled through the source term.

In this instance I found it best to solve the system iteratively using the TRI-DIAGONAL MATRIX ALGORITHM, a version of which is included in the code below.

Code:

import numpy as np
import matplotlib.pyplot as plt


#-----------------------------------------------------------


def tdma( a, b, c, d ):
#******************************************************************************
# Solve, using the TRI-DIAGONAL MATRIX ALGORITHM (TDMA)
#     a_i x_i-1 + b_i x_i + c_i x_i+1 = d_i,     i = 0, n-1
#
# Effectively, this is the n x n matrix equation, where
# a[i], b[i], c[i] are the non-zero diagonals of the matrix and d[i] is the rhs.
# a[0] and c[n-1] are not used.
#******************************************************************************
    n = len( d )
    p = np.zeros( n );   q = np.zeros( n )
    x = np.zeros( n )

    # Forward pass
    i = 0
    denominator = b[i]
    p[i] = -c[i] / denominator
    q[i] =  d[i] / denominator
    for i in range( 1, n ):
        denominator = b[i] + a[i] * p[i-1]
        if ( abs( denominator ) < 1.0e-10 ):
            print( "No solution" )
            return x
        p[i] =  -c[i]                   / denominator
        q[i] = ( d[i] - a[i] * q[i-1] ) / denominator

    # Backward pass
    x[n-1] = q[n-1]
    for i in range( n - 2, -1, -1 ):
        x[i] = p[i] * x[i+1] + q[i]

    return x


#-----------------------------------------------------------


N = 400                                # number of cells (1-N); added boundaries 0 and N+1
rmin, rmax = 0.0, 20.0                 # minimum and maximum r
dr = rmax / N                          # cell width
fL, fR, hL, hR = 0, 1, 0, 1            # boundary conditions on f and h at L (left) and R (right) boundaries
lamda = 0.5


# 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 nodes (except for boundaries)
r[0] = rmin;   r[N+1] = rmax                                     # boundary nodes on faces of cells
rv = np.linspace( rmin, rmax, N + 1 )                            # cell faces (index is the right side of the cell)


# Connection coefficients
awf = np.zeros( N + 2 )
aef = np.zeros( N + 2 )
awh = np.zeros( N + 2 )
aeh = np.zeros( N + 2 )
# Main cells
aef[1:N+1] = 1              / dr ** 2
aeh[1:N+1] = rv[1:N+1] ** 2 / dr ** 2
awf[1:N+1] = aef[0:N]
awh[1:N+1] = aeh[0:N]
# Boundaries
awf[1], awh[1], aef[N], aeh[N] = 0, 0, 0, 0     # zero out the matrix coefficients
awfL = 2             / dr ** 2                  # apply boundary conditions via the source term
aefR = 2             / dr ** 2
awhL = 2 * rmin ** 2 / dr ** 2
aehR = 2 * rmax ** 2 / dr ** 2


# Initial values
f = np.zeros( N + 2 )
h = np.zeros( N + 2 )
# Boundary conditions
f[0] = fL;   f[N+1] = fR;   h[0] = hL;   h[N+1] = hR


alpha = 0.5                  # Under-relaxation factor
niter = 0
s  = np.zeros( N + 2 )
sp = np.zeros( N + 2 )
for _ in range( 2000 ):
    niter += 1

    ftarget = f.copy()
    htarget = h.copy()

    ######################
    # Update f equation  #
    ######################
    # Source term (write as S + Sp.f, where Sp is negative)
    for i in range( 1, N + 1):                             # cells 1 - N
       s [i] = + h[i] ** 2 / 4
       sp[i] = - h[i] ** 2 / 4
       if ( 1 - f[i] ) * ( 1 - 2 * f[i] ) > 0:
          sp[i] += -2        * ( 1 - f[i] ) * ( 1 - 2 * f[i] ) / r[i] ** 2
       else:
          s [i] += -2 * f[i] * ( 1 - f[i] ) * ( 1 - 2 * f[i] ) / r[i] ** 2
    # Boundary conditions (applied through the source term)
    s[1] += awfL * fL;   sp[1] -= awfL
    s[N] += aefR * fR;   sp[N] -= aefR
    # Diagonal term (including implicit part of the source)
    ap = awf + aef - sp
    # Solve tri-diagonal system
    ftarget[1:N+1] = tdma( -awf[1:N+1], ap[1:N+1], -aef[1:N+1], s[1:N+1] )


    ######################
    # Update h equation  #
    ######################
    # Source term (write as S + Sp.h, where Sp is negative)
    for i in range( 1, N + 1):                             # cells 1 - N
       s [i] =                             lamda * r[i] ** 2 * h[i]
       sp[i] = -( 2 * ( 1 - f[i] ) ** 2 +  lamda * r[i] ** 2 * h[i] ** 2 )
    # Boundary conditions (applied through the source term)
    s[1] += awhL * hL;   sp[1] -= awhL
    s[N] += aehR * hR;   sp[N] -= aehR
    ap = awh + aeh - sp
    # Solve tri-diagonal system
    htarget[1:N+1] = tdma( -awh[1:N+1], ap[1:N+1], -aeh[1:N+1], s[1:N+1] )


    change = ( np.linalg.norm( ftarget - f ) + np.linalg.norm( htarget - h ) ) / N

    # Under-relax the update
    f = alpha * ftarget + ( 1 - alpha ) * f
    h = alpha * htarget + ( 1 - alpha ) * h

    # print( niter, change )
    if change < 1.0e-10: break

print( niter, " iterations " )


plt.plot( r, f, label='f' )
plt.plot( r, h, label='h' )
plt.xlim( ( 0, 10  ) )
plt.ylim( ( 0, 1.1 ) )
plt.grid()
plt.legend()
plt.show()

Output:

Compare Figure 11.6 in your linked document: