# How to calculate smoothing spline coefficients

I am attempting to calculate smoothing spline coefficients based on the description in Reinsch's 1967 paper, but I'm having some trouble. The first derivative is not continuous. Here are the equations I'm attempting to follow:

And this is the MATLAB code. There are some simplifications because all h_sub_i and weight values are equal to 1.

pts = [0      0
1      1
2      2
3      3
4      2
5      1
6      0];

n = length(pts)-1;
p = 10; % closeness of fit

Q = -2*eye(n-1) + diag(ones(1,n-2),1) + diag(ones(1,n-2),-1);
Q(n+1, end) = 0;
Q(n, end) = 1

T = 4/3*eye(n-1) + diag(1/3*ones(1,n-2),1) + diag(1/3*ones(1,n-2),-1)
Ap = transpose(Q)*Q + p*T

c = Ap^-1 * p * transpose(Q) * pts(:, 2);
a = pts(:, 2) - (1/p)*Q*c;
c = [0; c; 0];
d = (c(2:end) - c(1:end-1)) / 3;
b = (a(2:end) - a(1:end-1)) - c(1:end-1) - d;

N = 10000;
breaks = 0:6;
pp = mkpp(breaks, [d(1:n) c(1:n) b(1:n) a(1:n)]);
x = 0:0.01:6;
y = ppval(pp, x);

figure(1)
clf
hold on
grid on
plot(pts(:,1), pts(:,2), '--bo');
plot(x, y, 'r');


The following output is produced:

Q =

-2     1     0     0     0
1    -2     1     0     0
0     1    -2     1     0
0     0     1    -2     1
0     0     0     1    -2
0     0     0     0     1
0     0     0     0     0

T =

1.3333      0.33333            0            0            0
0.33333       1.3333      0.33333            0            0
0      0.33333       1.3333      0.33333            0
0            0      0.33333       1.3333      0.33333
0            0            0      0.33333       1.3333

Ap =

18.333     -0.66667            1            0            0
-0.66667       19.333     -0.66667            1            0
1     -0.66667       19.333     -0.66667            1
0            1     -0.66667       19.333     -0.66667
0            0            1     -0.66667       19.333


What is my error?