#!/usr/bin/env python3
import matplotlib.pyplot as plt

def p(x):
    return 2.0

def q(x):
    return 1.0 / x

def f(x):
    return -3.0

sigma1, sigma2 = 1.0, 0.5
gamma1, gamma2 = 0.0, -1.0
delta1, delta2 = 2.0, 1.0


def solve():
    N = 100
    l, r = 0.2, 0.5
    h = (r - l) / N

    def x(i):
        return l + i * h

    F = [-2.0 * h * h * f(x(i))         for i in range(N + 1)]
    A = [2.0 - p(x(i)) * h              for i in range(N + 1)]
    C = [2.0 * q(x(i)) * h * h - 4.0    for i in range(N + 1)]
    B = [2.0 + p(x(i)) * h              for i in range(N + 1)]

    p0 = gamma1 / (gamma1 - sigma1 * h)
    q0 = h * delta1 / (sigma1 * h - gamma1)
    pn = gamma2 / (gamma2 + sigma2 * h)
    qn = h * delta2 / (sigma2 * h + gamma2)

    F[1] -= A[1] * q0
    C[1] += A[1] * p0
    F[N - 1] -= B[N - 1] * qn
    C[N - 1] += B[N - 1] * pn

    alpha = [0, p0]
    beta = [0, q0]

    for i in range(1, N):
        alpha.append(-B[i] / (A[i] * alpha[i] + C[i]))
        beta.append((F[i] - A[i] * beta[i]) / (A[i] * alpha[i] + C[i]))

    y = [0] * (N + 1)
    y[N] = (qn + pn * beta[N]) / (1.0 - pn * alpha[N])

    for i in range(N - 1, -1, -1):
        y[i] = alpha[i + 1] * y[i + 1] + beta[i + 1]

    xl = [x(i) for i in range(N + 1)]
    plt.plot(xl[1:], y[1:])
    plt.show()

solve()