Question: 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

Answer: 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 r2. Collecting the derivative terms you will then get equations which look like spherically-symmetric diffusion equations:  Here, Sf and Sg 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 [ri-1/2,ri+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()