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()