How to solve this 2nd ODE with simple inequality constraints numerically? I have a 2nd ODE (derived from an elastic rod deflated naturally under its self-weight):
$$
y'' + K (x - 1) \cos(y) = 0 
$$
K is a constant coefficient. The variables range between $x \in [0, 1], y \in [-\pi / 2, 0]$
The BCs:
$$
y(0) = 0 \\
y'(1) = 0
$$
and a constraint $y' \le 0$ and $y'' \ge 0$ (from physical meaning).

How can I solve it numerically?
I have a weak background in the theory/classfication of ODEs, but I still tried to analysis and solve it as followed:
1. My own analysis on this problem
This is an second order nonlinear ODE. Its BCs are located in two sides, so it is a BVP.
A common practice is to split it into a 1st ODE system:
$$
\begin{cases}
y_2 = y_1' \\
y_2' + K (x - 1) \cos (y_1) = 0
\end{cases}
$$
Then I discrete all derivates by explicit euler method, and employe shooting method(where I guess the init BCs y'(0) and adjust this guess for a satisfying solution).
When K is big, I can solve it properly. But when K is small, my shooting method always failed. I have no idea how to analysis it anymore.
2. My questions

*

*Is my judgement about this problem correct? Is it an 2nd order ODE, Bounday value problems with not enough BCs?

*Why I can solve it properly in non-stiff case, and fail in stiff case?

*Some literatures point out that, auto07p package can help me on this problem. But I cannot understand its manual even a bit. Which book should I read to handle this problem?

*How can I know whether this problem have a solution?

*How to solve this problem numerically?

Sorry for my weak background in ODE.
Thanks a lot for your time!
3. solve_bvp and shooting failed when K > 200 and ignore the constraint
When K > 200, we need to set up a special init guess for solve_bvp, otherwise we will fail. Please check @Lutz 's answer below.
$

 A: In python you do
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt
from numpy import linspace, logspace, exp, cos

K = 5

ode = lambda x,y: [ y[1], K*(1-x)*cos(y[0]) ]
bc = lambda y0,y1: [ y0[0], y1[1] ]

x = linspace(0,1,5)
y = [0*x]*2
res = solve_bvp(ode, bc, x, y, tol=1e-5)

print(res.message)

x = linspace(0,1,150)
y = res.sol(x)
plt.plot(x,y[0])
plt.plot(res.x,res.y[0],'x')
plt.grid(); plt.show()

The solver uses a multi-shot approach where the segments are encoded using a collocation method. This produces a huge non-linear system that is solve via Newton using sparse linear algebra.

For large $K$ you get a nearly singular DE, so that a boundary layer approach becomes advisable for the initial guess. The outer solution is $y=-\frac{\pi}2$ (or any other root of the cosine). The inner solution at the left boundary with $Y(s)=y(K^{-1/2}s)$ follows the approximate DE
$$
Y''(s)=K^{-1}y''(K^{-1/2}s)\approx\cos(Y(s))
$$
This has a stable solution that approximately looks like $\frac\pi2(e^{-s}-1)$. So put that as initial guess
x = logspace(-6,0,21); x[0]=0;
y = [ pi/2*(exp(-x*K**0.5)-1),pi/2*K**0.5*exp(-x*K**0.5) ]

and I still get a solution for $K=500$,

$K=5000$ works too. It is also sufficient to simplify the initial guess to
x = logspace(-6,0,21); x[0]=0;
y = [0*x-pi/2, 0*x]

