q.py

import numpy as np
import scipy
from scipy.integrate import solve_bvp

P = 0.49

def rhs(x, y, p):
    λ1 = p[0]
    return np.vstack((y[1, :], -np.pi*np.cos(np.pi*y[0, :]) / (2 * λ1), y[1, :]**2))

def bc(ya, yb, p):
    return np.array([ya[0], yb[0] - 0.5, ya[2], yb[2] - P])

def main():
    x0 = np.linspace(0, 0.5, 32)
    y0 = np.vstack((x0, np.ones_like(x0), x0))
    sol = solve_bvp(rhs, bc, x0, y0, p=np.array([0.23]))
    print(sol)
    print(sol.success)
    return sol