% Function implementing natural spline (follow algorith 3.4)
function natural_spline
%nodes=[0.5, 1, 2, 3, 4, 8,10];
%nodes= [0:.2:10];
nodes=[0.1,0.2,0.3,0.4];
%fnodes=f(nodes);
fnodes=[-0.62049958, -0.28398668, 0.00660095, 0.24842440];
%%%%%%%%%%%%%%%%%%%%%%%
% Calculate Coefficients of the natural spline 
% given vaue of the function (fnodes) and the nodes
%%%%%%%%%%%%%%%%%%%%%%% 
n=size(nodes); n=n(2);               %determine number of the nodes
h = nodes(2:n)-nodes(1:n-1);     %create space for meshsteps : h(i)=nodes(i+1)-nodes(i);
a=fnodes;                      % calculate coefficients a_{i} (coinside with the value of the function at nodes)
c=zeros(1, n);  % setup space for coefficients c
b=zeros(1, n-1);  % setup space for coefficients b
d=zeros(1,n-1);   % setup space for coefficients d
%%%%%%%%%%%%%%%%%%%%%%
% calculate coefficients c_{i}, b_{i}, d_{i}
%%%%%%%%%%%%%%%%%%%%%%
alpha=(3./h(2:n-1)).*(a(3:n)-a(2:n-1))-(3./h(1:n-2)).*(a(2:n-1)-a(1:n-2)); %% calculate the right side of the system of equations for coefficients c_{i}
%% Solving a tridiagonal linear system by algorithm 6.7
l=zeros(1,n); l(1)=1; mu=zeros(1,n-1); z=zeros(1,n-1); % set up coefficients l, mu, z (step 3)
for i=2:1:n-1
  l(i)=2*(nodes(i+1)-nodes(i-1))-h(i-1)*mu(i-1);
  mu(i)=h(i)/l(i);
  z(i)=(alpha(i-1)-h(i-1)*z(i-1))/l(i);
end 
%% Step 5:
l(n)=1; c(n)=0; %z(n)=0
%% Step 6:
for j=n-1:-1:1
 c(j)=z(j)-mu(j)*c(j+1);
end
for j=1:1:n-1
 b(j)=(a(j+1)-a(j))/h(j)-h(j)*(c(j+1)+2*c(j))/3;
 d(j)=(c(j+1)-c(j))/(3*h(j));
end
spline(nodes,a,b,c,d,0.25)
spline_der(nodes,b,c,d,0.25)
%%%%%%%%%%%%%%%%%%%%%%%%%
% Plot the original function and the spline
%%%%%%%%%%%%%%%%%%%%%%%
figure(1);
subplot(1,1,1);
clf; % clear the plot
% Plot the exact function 
x=[0:.01:10];
%plot(x, f(x), '-r');
%hold on; % to keep the plot on 
% calculate the spline  
m=size(x);
m=m(2);
y=zeros(1,m);
for i=1:1:m
y(i)=spline(nodes,a,b,c,d,x(i));
end
plot(x,y-f(x),'-b');
hold off;


%%%%%%%%%%%%%%%%%%%%%%%%%
% at each interval [x_{j} ,x_{j+1}] spline is calculated by the formula
% a(j) + b(j)(x-x(j)) + c(j)(x-x(j))^2 + d(j)(x-x(j))^3
%  is the point where the spline needs to be evaluated
% %%%%%%%%%%%%%%%%%%%%%%
function y=spline(nodes,a,b,c,d,x);
y=0;
n=size(nodes); n=n(2);               %determine number of the nodes
for j=1:1:n-1
   if (x>=nodes(j)) & (x<=nodes(j+1)) 
      y=a(j) + b(j)*(x-nodes(j)) + c(j)*(x-nodes(j))^2 + d(j)*(x-nodes(j))^3;
   end
end  

function y=spline_der(nodes,b,c,d,x);
y=0;
n=size(nodes); n=n(2);               %determine number of the nodes
for j=1:1:n-1
   if (x>=nodes(j)) & (x<=nodes(j+1)) 
      y= b(j) + 2*c(j)*(x-nodes(j)) + 3*d(j)*(x-nodes(j))^2;
   end
end  

function y=f(x)
%y=log(x); 
y=sin(x-5);