Trying to obtain coefficient matrix of Bspline matrix then obtaining its quasi-diagonal matrix ?

I tried to obtain the coefficients of B-spline matrix using the following code :

function C_matrix = calculateBsplineCoefficientMatrix(s)
% Number of input samples
N = length(s);
% Construct Alpha Values (Used in the Recursive Calculation)
alpha_values = zeros(1, 10);
alpha_values(1) = -1/4;
for i = 2:10
    alpha_values(i) = -1 / (4 + alpha_values(i - 1));
end
alpha = alpha_values(end);
b_i = -alpha / (1 - alpha^2);
% Construct Forward Recursion Matrix (C_plus)
C_plus = zeros(N, N); % Matrix to store forward recursion values
C_plus(:, 1) = 1; % Initialize first column as 1
for k = 2:N
    C_plus(k, :) = alpha * C_plus(k-1, :);
    C_plus(k, k) = 1; % Set diagonal values to 1 for recursion steps
end
% Construct Backward Recursion Matrix (C_minus)
C_minus = zeros(N, N);
I_N = eye(N);  % Identity matrix of size NxN
C_minus(N, :) = b_i * (2 * C_plus(N, :) - I_N(N, :));
for k = N-1:-1:1
    C_minus(k, :) = alpha * (C_minus(k+1, :) - C_plus(k, :));
end
% Construct Final B-spline Coefficient Matrix
C_matrix = 6 * C_minus;
end

s here represents the signal or sequence basically this code is created from the paper FastO-line interpolation by Dooley in pdf document(1).

s = rand(1, 4);
C_spline_matrix = calculateBsplineCoefficientMatrix(s);
disp('B-Spline Coefficients Matrix Form:');
disp(C_spline_matrix);

The Solution is B-Spline Coefficients Matrix Form:

B-Spline Coefficients Matrix Form:

1.7327 -0.4665 0.1333 -0.0333

-0.4665 1.7410 -0.4974 0.1244

0.1333 -0.4974 1.8564 -0.4641

-0.0666 0.2487 -0.9282 1.7321

Then I tried to calculate its quasi diagonal matrix using the following code :

M = length(s); % Polynomial order or interpolation order
% Compute Transformation Matrices
T_mu = compute_TD(M); % Transformation matrix T_mu
T_z = calculate_Tz(M); % Transformation matrix T_z
% Ensure numerical stability
T_mu_invT = inv(transpose(T_mu)); % Explicit inverse transpose of T_mu
T_z_inv = inv(T_z);     % Explicit inverse of T_z
% Compute the quasi-diagonal matrix C_LCN
C_LCN = T_mu_invT .*C_spline_matrix .* T_z_inv;
% Display the final quasi-diagonal matrix
disp('Computed C_LCN Matrix:');
disp(C_LCN);

The solution of quasi diagonal matrix

Computed C_LCN Matrix:

1.8384 0 0 0

-0.5937 -2.4397 0 0

0.1728 1.1772 0.8176 0

-0.0168 -0.2263 -0.3135 -0.0388

The relevant functions are given below:

function Td1 = compute_Td1(M)
% Create a matrix of binomial coefficients
[I, J] = ndgrid(1:M, 1:M);
binomials = zeros(M, M);
for ii = 1:M
    for jj = 1:M
        if jj <= ii
            binomials(ii, jj) = nchoosek(ii-1, jj-1);
        end
    end
end
% Compute powers and combine with binomial coefficients
powers = ((- (M - 1) / 2) .^ (I - J)) .* (J <= I);
Td1 = binomials .* powers;
% Debugging: Display Td1 matrix
disp('Td1 Matrix:');
disp(Td1);
end
function Td2 = compute_Td2(M)
% Compute the Td2 matrix
Td2 = eye(M+1);
for n = 2:M+1
    for k = 2:n-1
        Td2(n, k) = Td2(n-1, k-1) - (n-2) * Td2(n-1, k);
    end
end
Td2 = Td2(2:M+1, 2:M+1);
% Debugging: Display Td2 matrix
disp('Td2 Matrix:');
disp(Td2);
end
function TD = compute_TD(M)
Td1 = compute_Td1(M);
Td2 = compute_Td2(M);
% Compute TD
TD = Td1 * Td2;
% Force symmetry in TD
TD = (TD + TD') / 2;
% Regularization for numerical stability
TD = TD + 1e-8 * eye(size(TD));
% Debugging: Display TD matrix
disp('Symmetrized TD Matrix:');
disp(TD);
end
function Tz = calculate_Tz(M)
% Generate the Tz transformation matrix for Newton interpolation
Tz = zeros(M, M);
for i = 1:M
    for j = 1:i
        Tz(i, j) = nchoosek(i-1, j-1) * (-1)^(j+1);
    end
end
% Regularization for numerical stability
Tz = Tz + 1e-8 * eye(size(Tz));
% Debugging: Display Tz matrix
disp('Tz Matrix:');
disp(Tz);
end

The problem is the it does not match the values as solved in the paper Electronic letters...

The solutions provided are :

Please could any body help me correcting this code to obtain the solution ?