How to compute the following Crank-Nicolson type scheme?

Question

Context: I am attempting to replicate the results of Crank-Nicolson Type Method for Burgers Equation .

From my understanding of Crank-Nicolson schemes, one can set up a tri-diagonal matrix and "conveniently" solve the system using the Thomas algorithm.

However, the discretized equation (4) is as follows.

$$ (-r+pu^{n}_j)u^{n+1}_{j+1}+(2r+1+pu^{n}_{j+1}-pu^{n}_{j-1})u^{n+1}_{j} -(r+pu^{n}_{j})u^{n+1}_{j-1}=ru^{n}_{j+1}+(1-2r)u^{n}_{j}+ru^{n}_{j-1} $$

How would I be able to set up the tri-diaganol matrix here? The equation can't seemed to be simplified further to have all the $j+1$ terms on one side while the $j$ terms on the other. Is there another method for me to compute this scheme?

Is this similar to Crank-Nicolson for quadratic PDE where I'm stuck with a non-linear equation?

I apologize if this seems trivial.

---

Update

(I'm not sure if it remains appropriate for me to continue asking here, to start a new question or perhaps better to ask in stack overflow, please let me know. The question remains to compute the Crank-Nicolson scheme so I'll continue asking here first.)

As mentioned in the comment, it should be $n$ to $n+1$, hence possible to set this as a tri-diagonal matrix.

I have now attempted to compute the scheme but it appears I have done something wrong when comparing it with known solutions. I suspect that I have not used the boundary and initial conditions correctly, or not followed the indices properly.

As included within the article, for $t>0$ and $x \in (0,1)$ the boundary conditions are:

$$u(x=0,t)=u(x=1,t)=0 $$

and initial condition for $0<x<1$:

$$u(x,0)=sin(\pi x)$$.

For your quick reference, in my code (in R) I have implemented the coefficients on the LHS of discretized equation(4) as A,B and C. D consists of the entire RHS. $$A_j=(-r+pu^{n}_j)$$ In $B_j$, I suspect something was not considered correctly when computing it, $$B_j=(2r+1+pu^{n}_{j+1}-pu^{n}_{j-1})$$
$$C_j=-(r+pu^{n}_{j})$$ and $$D_j=ru^{n}_{j+1}+(1-2r)u^{n}_{j}+ru^{n}_{j-1}$$

---

Update 2

I have found some mistakes in the implementation where I did not consider the initial and boundary conditions correctly. I have updated the code accordingly. However, I am still unable to obtain the correct results.

# Some parameters as given in the article
N = 1/0.025
input_dx = 0.025
k=10
input_dt=0.0001
p=input_dt/(4*input_dx)
r=(k*input_dt/(2*(input_dx^2)))
time_inter=81

# 0 to N. So If N=40, expect 41 values.
# IC for 0< j <N . In R, # equivalently: 1< j <N+1. 
# Therefore 2:N
# Spatial discretization
x_discrete=numeric()
x_discrete[1]=0
for( j in 2:(N)){
    x_discrete[j] = x_discrete[1] + (j-1)*input_dx
}
x_discrete[N+1]=1


# Create empty matrix to store solutions
u_sol = matrix(0, (N+1), (time_inter+1))

# BC u(0,t>0) = u(1, t>0) = 0 (Redundant since matrix is made up of all zeroes, but included)
u_sol[1,2:time_inter]=0
u_sol[N+1,2:time_inter]=0

# IC
u=numeric()
for( j in 1:length(x_discrete)){
    u[j] = sin(pi*x_discrete[j])
    }

u_sol[,1]=u

# For each time-step
for(n in 1:time_inter){
    # A,B,C, D set up from:
    # https://www.sfu.ca/~rjones/bus864/notes/notes2.pdf        
    # RHS of equation
    D=numeric()
    for(j in 2:N){
        if(n==1){
            # When n=0, use IC
            D[j-1] = r*u_sol[(j+1), n] + (1-2*r)*u_sol[(j), n] + r*u_sol[(j-1), n]
        }else if(n>1){
            if(j==2){
                # When j=1, implement BC as j-1=0
                D[j-1] = r*u_sol[(j+1), n] + (1-2*r)*u_sol[(j), n] + 0 # Adjust according to BC
            }else if(j==N){
                # When j=N-1, implement BC as j+1=N 
                D[j-1] = 0 + (1-2*r)*u_sol[(j), n] + r*u_sol[(j-1), n] # Adjust according to BC
            }else if((j>2)&(j<(N))){
                # When 1<j<N-1 (Note: N-1 in comment means N-2 index in R)
                D[j-1] = r*u_sol[(j+1), n] + (1-2*r)*u_sol[(j), n] + r*u_sol[(j-1), n]
            }
        }
    }
    # LHS of equation
    A=numeric() # Coefficient of u^{n+1}_{j+1}
    B=numeric() # Coefficient of u^{n+1}_{j}
    C=numeric() # Coefficient of u^{n+1}_{j-1}
    for(j in 2:(N)){
        if(n==1){
            # When n=0, all coefficients for all x values are calculated via IC
            A[j-1] = -r + p*u_sol[(j), n]
            B[j-1] = 2*r + 1 + p*u_sol[(j+1), n] - p*u_sol[(j-1), n] # As seen in article; equation (4)
            C[j-1] = -(r + p*u_sol[(j), n])
        }else if(n>1){
            # When n>0
            if(j==N){
                # If j=N-1, j+1=N. Hence, use BC
                A[j-1] = -r + p*u_sol[(j), n]
                B[j-1] = 2*r + 1 + p*0             - p*u_sol[(j-1), n] # Adjust according to BC
                C[j-1] = -(r + p*u_sol[(j), n])         
            }else if(j==2){
                # If j=1, j-1=0. Hence, use BC
                A[j-1] = -r + p*u_sol[(j), n]
                B[j-1] = 2*r + 1 + p*u_sol[(j+1), n] - p*0             # Adjust according to BC
                C[j-1] = -(r + p*u_sol[(j), n])                     
            }else{      
                A[j-1] = -r + p*u_sol[(j), n]
                B[j-1] = 2*r + 1 + p*u_sol[(j+1), n] - p*u_sol[(j-1), n]
                C[j-1] = -(r + p*u_sol[(j), n])
            }
        }
    }

    # Set up tridiagonal-matrix
    # Method taken from:
    # https://stackoverflow.com/questions/49536400/create-a-tridiagonal-matrix
    m = diag(B)
    m[row(m) - col(m) == 1] = A[-1]
    m[row(m) - col(m) == -1] = C[-length(C)]

    # Solve for solutions
    sol=MASS::ginv(m)%*%t(t(D))

    # Update solution matrix
    u_sol[2:N, (n+1)]=sol
}


# To check, at 80 for k=10, x=0.1, dt=0.0001
# Exact solution=0.11461 compare with
u_sol[5, 81]