How to solve a non linear ODE with Newton's method? Im trying to solve this ODE using the Newton method (for non-linear equations using Jacobian Matrix):
$$
\frac{u''(x)}{(1+u'(x)^2)^{3/2}}=\frac TK u(x)+\frac{wx(x-L)}{2K}
$$
original ode image
$u(0) = u(L) = 0$ (let's take $L=9)$
$T=500$
$K= 5 \times 10^9$
$w = 100$
$u''$ and $u'(x)$ are defined as follows:

I don't know what is f1 and what is f1. How could this be implemented in Python? Anyone can help here please?
 A: You get an automatic Newton solver by using the BVP solver scipy.integrate.solve_bvp. You program your equation as usual
def eqn(x,u): return [ u[1], (1+u[1]**2)**1.5*(T*u[0]+0.5*w*x*(x-L))/K ]
def bc(u0,uL): return [ u0[0], uL[0] ]

x = np.linspace(0,L,21);
a = 0.1
u = [ a*x*(L-x), a*(L-2*x) ]
res = solve_bvp(eqn, bc, x, u, tol = 1e-12); # the solution is small, scale 1e-6, thus reduce from the default tol=1e-3
print res.message
if res.success:
    x = np.linspace(0,L,700)
    u = res.sol(x)
    plt.plot(x,u[0])

resulting in the plot

A: Given a second order boundary value problem $y''(x)=f(x,y(x),y'(x))$, $y(a)=y_a$, $y(b)=y_b$, you can linearize it at a given approximate solution $y$ as
\begin{align}
y''(x)+Δy''(x) &= f(x,y(x)+Δy(x),y'(x)+Δy'(x))
\\
&\approx f(x,y(x),y'(x))+\partial_yf(x,y(x),y'(x))Δy(x)+\partial_{y'}f(x,y(x),y'(x))Δy'(x)
\end{align}
Now insert the difference quotients for $y$ to produce the coefficients and use the difference quotients for $Δy$ to establish the usual linear system for
$$
Δy''(x) + a(x)Δy'(x)+b(x)Δy(x)=c(x)
$$
leading to
$$
(1+\tfrac12ha_k)Δy_{k+1}-(2-h^2b_k)Δy_k+(1-\tfrac12ha_k)Δy_{k-1}=c_k
\\~\\
a_k=-\partial_{y'}f(x,y_k,\frac{y_{k+1}-y_{k-1}}{2h})
\\
b_k=-\partial_{y}f(x,y_k,\frac{y_{k+1}-y_{k-1}}{2h})
\\
c_k=f(x,y_k,\frac{y_{k+1}-y_{k-1}}{2h}).
$$

With the given problem we have
$$
f(x,y,y')=\frac{Ty-\tfrac12wx(L-x)}K(1+y'^2)^{3/2}
$$
so that
$$
\partial_{y}f=\frac{T}K(1+y'^2)^{3/2},\\
\partial_{y'}f=\frac{Ty-\tfrac12wx(L-x)}K\cdot 3(1+y'^2)^{1/2}y'
$$
which now needs to be assembled into a code establishing the linear system for the $Δy_k$, which then is iteratively solved and used to update $y_k$.
