Natural cubic spline as piecewise linear + correction term

naturalcubic.m

function naturalcubic
a = 1995;
b = 2000;
y = [30 25 20 28 15 31];
n = length(y);
h = (b-a)/(n-1);
rhs = (-6/h^2)*[0; diff(y',2); 0];
A = diag(4*ones(1,n))+diag(ones(1,n-1),1)+diag(ones(1,n-1),-1);
A(1, 1:2) = [1 0];
A(n, n-1:n) = [0 1];
plotspline(a, b, y, A\rhs, 'r');
hold on
plot(a:h:b, y, 'r*'); 
end 

function plotspline(a,b,y,z,color)
h = (b-a)/(length(y)-1);  
allx = [];
spline = [];
for i = 1:length(y)-1  
   xL = a+h*(i-1);
   xR = a+h*i;
   x = linspace(xL, xR);
   linear = y(i)*(xR-x)/h + y(i+1)*(x-xL)/h;
   correction = ((z(i+1)+2*z(i))*(xR-x)+(2*z(i+1)+z(i))*(x-xL)).*(xR-x).*(x-xL)/(6*h);
   allx = [allx, x];   
   spline = [spline, linear+correction]; 
   plot(x, linear, 'g');
   plot(x, correction, 'm');
end
plot(allx, spline, color);
end