# code to solve, and then plot, the solution of the nonlinear BVP
#        y'' = f(x, y, y')     for xL < x < xR
# where
#          y(xL) = yL ,  y(xR) = yR

# an explanation of how to use this code is given at the end of this file

import numpy as np
import matplotlib.pyplot as plt


def nbvp():
    # set boundary conditions
    xL = 0
    yL = 3
    xR = 1
    yR = 5

    # parameters for calculation
    nx = 200
    error = 0.000001

    # start off with a linear solution
    x = np.linspace(xL, xR, nx + 2)
    y = np.zeros(nx + 2)
    for ix in range(nx + 2):
        y[ix] = yL + (yR - yL)*x[ix]

    # plt.clf()
    # fig = plt.gcf()
    # fig.set_size_inches(5.73, 1.99)

    dx = x[1] - x[0]
    dxx = dx*dx
    err = 1

    counter = 0
    while err > error:

        # calculate the coefficients of finite difference equation
        a = np.zeros(nx)
        c = np.zeros(nx)
        v = np.zeros(nx)
        u = np.zeros(nx)

        for j in range(1, nx + 1):
            jj = j - 1
            z = (y[j + 1] - y[j - 1])/(2*dx)
            a[jj] = 2 + dxx*fy(x[j], y[j], z)
            c[jj] = -1 - 0.5*dx*fz(x[j], y[j], z)
            v[jj] = (
                -2*y[j] + y[j + 1] + y[j - 1]
                - dxx*f(x[j], y[j], z)
            )

        # Newton iteration
        v[0] = v[0]/a[0]
        u[0] = -(2 + c[0])/a[0]

        for j in range(1, nx):
            xl = a[j] - c[j]*u[j - 1]
            v[j] = (v[j] - c[j]*v[j - 1])/xl
            u[j] = -(2 + c[j])/xl

        vv = v[nx - 1]
        y[nx] = y[nx] + vv
        err = abs(vv)

        for jj in range(nx - 1, 0, -1):
            vv = v[jj - 1] - u[jj - 1]*vv
            err = max(err, abs(vv))
            y[jj] = y[jj] + vv

        counter = counter + 1

    newton_iterations = counter
    print(f"newton_iterations = {newton_iterations}")

    # plot computed solution
    plt.plot(x, y)
    plt.grid(True)
    # plt.axis([-0.02, 1.02, -2, 2])
    plt.xlabel("x-axis", fontsize=14, fontweight="bold")
    plt.ylabel("Solution", fontsize=14, fontweight="bold")
    plt.tick_params(axis="both", labelsize=14)
    plt.box(True)
    plt.show()

    return x, y


#  right hand side of differential equation
def f(x, y, z):
    g = -1 + 20*y - np.tanh(z)
    return g


#  fy(x,y,z) = diff(f,y)
def fy(x, y, z):
    g = 20
    return g


#  fz(x,y,z) = diff(f,z)
def fz(x, y, z):
    g = -1/np.cosh(z)**2
    return g


if __name__ == "__main__":
    nbvp()


#  This code is run by entering the command
#
#      python nbvp.py
#
#  in a terminal or command window (making sure that this file is in the
#  current folder or that its folder is on the Python path). The user needs
#  to first enter certain information about the problem into this file before
#  running it.
#
#  The code requires NumPy and Matplotlib. If needed, install them with
#
#      python -m pip install numpy matplotlib
#
#  Needed Input From User
#
#  Boundary conditions (BCs): xL, yL, xR and yR in nbvp()
#  Number of grid points: nx in nbvp()
#
#  f(x,y,z):  this is the right hand side of the differential equation
#             with z acting as the place holder for y'
#  fy(x,y,z): this is the derivative of f with respect to y
#  fz(x,y,z): this is the derivative of f with respect to z
#
#  Example:   y'' = -1 + 20y - tanh(y') for 0 < x < 1
#             with y(0) = 3 and y(1) = 5
#    nx = 200
#    xL = 0, yL = 3, xR = 1, yR = 5
#    f(x,y,z) = -1 + 20*y - tanh(z)
#    fy(x,y,z) = 20
#    fz(x,y,z) = -sech(z)^2
#
#  What to do if it doesn't work (aside from checking on input errors)
#  1.  The program uses a linear function that satisfies the BCs as a start-off
#      guess. You might try a quadratic function like
#
#      y = (yL*(x - xR)**2 + yR*(x - xL)**2)/(xR - xL)**2
#
#  2.  If the problem has a boundary or interior layer you should consider
#      increasing nx (don't be shy about making nx very large).
#  3.  The error value in nbvp() is used by Newton's method. It is very
#      doubtful that you will need to change this, but certainly you can if
#      you feel that it will help.
#
#  Background
#  This is a simple code that solves nonlinear 2nd order BVPs. It does not
#  use adaptive methods, but instead uses 2nd order centered differences on
#  a uniform grid, with Newton's method to solve the resulting nonlinear system.
#  The specific algorithm is described in the book Intro to Numerical Methods
#  for Differential Equations (Holmes, 2007) on pages 59-62. It is very robust,
#  and was used to solve all of the nonlinear BVPs in Intro to Perturbation
#  Methods (Holmes, 1995). It is also relatively fast. The execution times
#  will depend on the Python installation, computer, and value of nx.
#
#  This Python version was translated from code tested on MATLAB R2009b.
#  Original MATLAB version: April 7, 2010