function err = CSpl(x,y,k,fy,yp)
%
%INPUT
%   x = vector of x data
%   y = vector of y=f(x) data
%   k = number of points to add to initial interpolation
%   fy = function value of interpolated points
%
%OUTPUT
%   err = sup err
%
%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%
%   Initializations
%
[m,n] = size(x);
h=diff(x);
M=diag(ones(n,1))+diag(4*ones(n-1,1),1)+diag(ones(n-2,1),2);
M = M(1:n-2,2:n-1);
s=6./(h(1:(end-1)).^2).*diff(diff(y)); s=s';
M=rref([M,s]);
Msol1=M(:,n-1);
Msol=(zeros(size(x)))';
Msol(2:n-1) = Msol1(1:n-2);
Msol(1) = 0;
Msol(n) = 0;  
a = zeros(1, n-1);
b = zeros(1, n-1);
c = zeros(1, n-1);
d = zeros(1, n-1);
%
%	Calculate the coefficients
%
for i = 1:(n-1)
   a(1,i) = (Msol(i+1)-Msol(i))/(6*h(1,i));
   b(1,i) = (Msol(i))/2;
   c(1,i) = (y(1,i+1)-y(1,i))/h(1,i)-(Msol(i+1)+2*Msol(i))*h(1,i)/6;
   d(1,i) = y(1,i);
end
%
%	Generate domain for Spline to be graphed 
%
Xi = zeros(n-1,k);
   for iter = 1:(n-1)
      domain = linspace(x(1,iter),x(1,iter+1),k);
      Xi(iter,:)=domain(1,:);
   end
%   
%	Plot functions and calculate error
%
plot(x,y,'ro');
hold on
Spl = zeros(k,n-1);
Fmat =zeros(k,n-1);
inferr = zeros(1,n-1);
%
%   Load function data in a matrix for computing error
%
for j = 1:n-1
    for i=1:k
        Fmat(i,j) = fy(i + k*(j-1));
    end
end
%
%   Calculate spline data and compute error
%
   for i = 1:(n-1)
      Smat = a(1,i)*(Xi(i,:)-x(i)).^3+ b(1,i)*(Xi(i,:)-x(i)).^2+c(1,i)*(Xi(i,:)-x(i))+ d(1,i);
      plot(Xi(i,:),Smat,'-','color', 'b')
      inferr(i) = max(abs(Fmat(:,i)' - Smat)/n);
      Spl(:,i) = Smat';
   end
err = max(inferr);
plot(yp,fy,'--','color','g')
legend('Data points','Spline approximation','Actual values')
title(err)
xlabel('l_{\infty} error');
hold off
end