Crank-Nicolson fix

I am trying to implement the Crank-Nicolson method, which is described in Burden and Faires 10th edition as below:

Image containing Crank-Nicolson Pseudo Algorithm

Observation: I was avoiding using images in the problem description, but could not think of a better way to present the steps.

I tried fixing the code in the related question Crank-Nicolson Method, but was unable to. However, I further commented the code and did find some mistakes, so I think I am closer to the correct implementation. The code is now as follows:

import numpy as np
np.set_printoptions(formatter={'float': '{: 0.3f}'.format})

def f(x):
    return np.sin(np.pi * x * 0.1)

def zero_matrix(i, j, start_from=1):
    matrix = [[-1] * (j+start_from)]
    for x in range(i):
        line = [-1] * start_from + [0] * j
        matrix.append(line)
    return matrix

def zero_vector(i, start_from=1):
   return [-1] * start_from + [0] * i

def algo_12_3(l, T, alpha, m, N):
     """Crank-Nicolson Method.

Args:
    l: endpoint
    T: maximum time
    alpha: constant
    m:
    N:
Return:
    approximations wi,j to u(xi, tj) for each i=i,...,m-1 and j=1,...,N.
"""
w = zero_matrix(m, N)
z = zero_vector(m-1)
ll = zero_vector(m-1)
u = zero_vector(m-1)

m = 10
k = 0.01
Lambda = 1
h = 0.1

#Step 2
for i in range(1, m):
    w[i][0] = f(i * h)
    
#Step 3
ll[1] = 1 + Lambda
u[1] = -Lambda/(2*ll[1])

for i in range(2, m-1):
    ll[i] = 1 + Lambda + (Lambda * u[i-1])/2
    u[i] = -Lambda/(2 * ll[i])
    ll[m-1] = 1 + Lambda + (Lambda * u[m-2])/2
    
#Step 6
for j in range(1, N + 1):
    
    #Step 7
    t = j * k
    
    #Step 8
    z[1] = ((1 - Lambda) * w[1][j - 1] + (Lambda / 2) * w[2][j-1]) / ll[1]
    
    for i in range(2, m):
        z[i] = ((1 - Lambda) * w[i][j-1] + (Lambda / 2) * (w[i+1][j-1] + w[i-1][j-1] + z[i-1])) / ll[i]
        
    #Step 9
    w[m-1][j] = z[m-1]
    
    #Step 10
    for i in range(m - 2, 1, -1):
        w[i][j] = z[i] - u[i] * w[i+1][j - 1]
        
    #Step 11
    print('t = {:.2f}: '.format(t))
    
    for i in range(1, m + 1):
        print('(xi = {:.2f}, wi,j = {:.8f})'.format(i * h, w[i][j]))

algo_12_3(1, 10, 1, 100, 50)

The last line of the above piece of code outputs:

Output

Which does not match the table found in the book:

Book output

I am also using the same parameters of the example. (Example 3, page 754, Chapter 12)

Can anyone please point some errors in the implementation?
Thanks in advance.

Source: Python-3x Questions

LEAVE A COMMENT