What causes the difference in y-axis plots in R and Python?

Question

Python

I rewrote the program from R to Python, but for some reason the y (x) values ​​​​on Python jump too much depending on N, with n = 100 the order of ordinates is dozens, with 1000 it is the order of hundreds, on R as it was within + - from 0 to 1 , and it remains, it turns out some kind of instability of the solution in python. Can someone explain what the issue is related to?
Here on R:

require(ggplot2)

ua = 0.8
ub = 1

A = function(x){
  return( log(3 + sqrt(1 + x) ))
}
B = function(x){
  return( exp(-x / 2) )
}
C = function(x){
  return ( sqrt((3 - x) / 2) )
}
F = function(x){
  return ( 2 - abs(1 - x) )
}

run = function(N){
  h = 2.4/(N + 1)
  
  an = function(x){
    return( -(2 * A(x) + B(x) * h) )
  }
  bn = function(x){
    return( 4 * A(x) + 2 * C(x) *h^2 )
  }
  cn = function(x){
    return( -(2 * A(x) - B(x) * h) )
  }
  fn = function(x){
    return( 2 * F(x) * h^2 )
  }
  
  x = seq(from = 0.1, to = 2.5, length.out = N+2)
  u = vector("numeric", N+2)
  a = sapply(x, an)
  b = sapply(x, bn)
  c = sapply(x, cn)
  f = sapply(x, fn)
  
  
  a[1] = 0
  b[1] = 1
  c[1] = 0
  f[1] = 0.8
  u[1] = 0.8
  
  a[N+2] = 0
  b[N+2] = 1
  c[N+2] = 0
  f[N+2] = 1
  u[N+2] = 1
 
  
  alpha = vector("numeric", N+2)
  beta  = vector("numeric", N+2)
  
  alpha[1] = -c[1]/b[1]
  beta[1]  = (f[1])/b[1]
  
  for(n in 2:(N + 1)){
    alpha[n]= -c[n] / (a[n] * alpha[n-1] + b[n])
    beta[n] = (f[n] - a[n] * beta[n-1]) / (a[n] * alpha[n-1] + b[n])
  }
  
  for(n in (N+1):2){
    u[n] = alpha[n] * u[n+1] + beta[n]
  }
  
  df = data.frame(x = x, y = u)
  g = ggplot(df, aes(x)) + geom_line(aes(y=y), col = "red") 
  print(g)
  print(u[N+2])
}
run(100)

Here it is in python:

import math as m
import matplotlib.pyplot as plt

def run(n):
    h = 2.4 / (n+1)

    a = list(range(n)); b = list(range(n)); c = list(range(n)); f = list(range(n))

    # начальные условия
    a[0] = 0; b[0] = 1; c[0] = 0 ; f[0] = 0.8
    a[-1] = 0; b[-1] = 1; c[-1] = 0; f[-1] = 1

    alpha = list(range(n)); beta = list(range(n))
    alpha[0] = (-1)*(c[0])/(b[0])
    beta[0] = (f[0])/(b[0])

    x = list(range(n))
    x[0] = 0.1

    for i in range(1,n):
        x[i] = x[0] + i*h

    A = list(range(n)); B = list(range(n)); C = list(range(n)); F = list(range(n))

    # функци A(x), B(x), C(x), f(x).
    for i in range(n):
        A[i] = m.log(3 + m.sqrt(1+x[i]))
        B[i] = m.exp(-x[i]/2)
        C[i] = m.sqrt((3-x[i])/2)
        F[i] = 2 - m.fabs(1-x[i])

    # коэф-ты. 
    for i in range(1,n):
        a[i] = (-1) * (2 * A[i] + h * B[i])
        b[i] = (4 * A[i] + 2 * pow(h, 2) * C[i])
        c[i] = (-1) * (2 * A[i] - h * B[i])
        f[i] = 2*pow(h, 2) * f[i]

        alpha[i] = (-1)*(c[i])/(a[i]*alpha[i-1]+b[i])
        beta[i] = (f[i] - a[i]*beta[i-1])/ (a[i]*alpha[i-1]+b[i])

    # Решение u_(N+1)
    y = list(range(n))
    y[0] = 0.8
    y[-1] = 1 # кладу последний n-й элемент равный 1
    for i in reversed(range(1, n-1)):
        y[i] = alpha[i]*y[i+1]+beta[i]

    plt.plot(x, y, color='red', marker='o', linestyle='solid')
    plt.show()

run(100)