Реализация метода прогонки в решении краевой задачи

У меня есть краевая задача

y'' + p(x) * y' + q(x) * y = f(x),
y(0) = y0,
y(L) = yN.

Надо решить её методом конечных разностей. Иными словами,

1. Ввести на заданном отрезке равномерную сетку,

2. Заменить во внутренних точках разбиения отрезка приближённо производные,

3. Подставить производные в исходное дифференциальное уравнение, получив систему линейных алгебраических уравнений с трёхдиагональной матрицей,

4. Решить получившуюся систему линейных алгебраических уравнений относительно функции y(x), учитывая начальные краевые условия

   y(0) = y0, y(L) = yN.

Но у меня есть проблема. Я добиваюсь выполнения условия

y(0) = y0,

но не условия

y(L) = yN.

Прикладываю код своей собственной программы на языке программирования Python:

 import sympy as sp
 import numpy as np
 import matplotlib.pyplot as plt

 x = sp.Symbol("x")

 y = sp.Function("y(x)")
 precise_decision = sp.Function("u(x)")

 p = sp.Function("p(x)")
 q = sp.Function("q(x)")
 f = sp.Function("f(x)")

 x0 = 0.0

 def launch(p, q, f, xN, y0, yN, N, precise_decision):
     if xN is np.nan or xN is np.inf or xN <= x0:
         print("Некорректно введён правый конец отрезка краевой задачи")
         return
     if y0 is np.nan or y0 is np.inf:
         print("Некорректно введено первое начальное условие краевой задачи")
         return
     if yN is np.nan or yN is np.inf:
         print("Некорректно введено второе начальное условие краевой задачи")
         return
     if N is np.nan or N is np.inf or N <= 0:
         print("Некорректно введено число отрезков, на которые разбивается отрезок [{0}; {1}]".format(x0, xN))
         return

     display_input_data(p, q, f, xN, y0, yN)

     step = (xN - x0) / N
     x_values = np.arange(x0, xN + step, step)
     y_values = tridiagonal_matrix_algorithm(p, q, f, step, N, x_values, y0, yN)
     paint_graphic(x_values, y_values, precise_decision)

     return

 def display_input_data(p, q, f, xN, y0, yN):
     sp.pprint(p)
     sp.pprint(q)
     sp.pprint(f)
     print("Заданный отрезок:    [{0}; {1}]".format(x0, xN))
     print("Начальные условия:   y(0) = {0}, y(l) = {1}".format(y0, yN))
     print()

     return

 def tridiagonal_matrix_algorithm(p, q, f, delta, N, x_values, y0, yN):
     alphas = []
     alpha = 0.0

     betas = []
     beta = 0.0
     for i in range(1, N + 1):
         x_curr = x_values[i]
         A = (1.0 / np.power(delta, 2)) - p.subs(x, x_curr).evalf() / (2.0 * delta)
         C = (1.0 / np.power(delta, 2)) + p.subs(x, x_curr).evalf() / (2.0 * delta)
         B = -1.0 * (A + C) + q.subs(x, x_curr).evalf()
         F = f.subs(x, x_curr).evalf()

         if i == 1:
             alpha = -1.0 * C / B
             beta = (F - A * y0) / B
         elif i < N:
             denominator = A * alpha + B
             alpha = -1.0 * C / denominator
             beta = (F - A * beta) / denominator
         else:
             alpha = 0.0
             beta = (F - B * yN) / A

         alphas.append(alpha)
         betas.append(beta)

     y_value = yN
     y_values = np.zeros((N + 1))
     y_values[N] = y_value
     for i in range(N - 1, -1, -1):
         y_value = alphas[i] * y_value + betas[i]
         y_values[i] = y_value
     y_values[0] = y0

     return y_values

 def paint_graphic(x_values, y_values, precise_decision):
     plt.figure(figsize = (12, 12))

     plt.plot(x_values, y_values, marker = "o", markersize = 6, color = "black", label = "приближённое решение y(x)", linewidth = 1.0, linestyle = "--")

     h = sp.lambdify(x, precise_decision, "numpy")
     mapping = np.linspace(x_values[0], x_values[len(x_values) - 1], 100000)
     plt.plot(mapping, h(mapping), label = "точное решение u(x)", linewidth = 2, linestyle = "-", color = "red")

     plt.xlabel("x", color = "blue", fontsize = 20)
     plt.ylabel("y(x)", color = "blue", fontsize = 20)

     plt.legend()

     plt.title("Графики приближённого и точного решений краевой задачи", color = "blue", fontsize = 20)

     plt.grid(linewidth = 1, linestyle = "--")
     plt.xticks(fontsize = 10)
     plt.yticks(fontsize = 10)

     plt.savefig("graphic.png", format = "png")

     plt.show()

     return

`

В качестве примеров были предложены краевые задачи:

p = 0.0 * x
q = 4.0 + 0.0 * x
f = 5.0 * sp.exp(x)
y(0.0) = 1
y(np.pi / 4) = np.exp(np.pi / 4) + 1.0
precise_decision = sp.sin(2.0 * x) + sp.exp(x)

и

p = sp.tan(x)
q = sp.Pow(sp.cos(x), 2)
f = 0.0 + 0.0 * x
y(0.0) = 1
y(8.0) = 10.0
precise_decision = sp.cos(sp.sin(x)) + ((10.0 - sp.cos(sp.sin(8))) / sp.sin(sp.sin(8))) * sp.sin(sp.sin(x))