function y_eval = ncs(x_nodes,y_nodes,x_eval)

m=length(x_nodes);
n=m-2;
a=zeros(1,n); b=zeros(1,n); c=zeros(1,n); f=zeros(1,n);
for i=1:n
    b(i)=(x_nodes(i+2)-x_nodes(i))/3;
end
for i=2:n
    a(i)=(x_nodes(i+1)-x_nodes(i))/6;
end
for i=1:n-1
    c(i)=(x_nodes(i+2)-x_nodes(i+1))/6;
end
for i=1:n
    f(i)=(y_nodes(i+2)-y_nodes(i+1))/(x_nodes(i+2)-x_nodes(i+1)) ...
            -(y_nodes(i+1)-y_nodes(i))/(x_nodes(i+1)-x_nodes(i));
end

[soln, alpha, beta, ier] = tridiag(a,b,c,f,n,0);

M=zeros(m);
M(2:m-1)=soln;

len_x=length(x_eval);
for i=1:len_x
    x=x_eval(i);
    for j=2:m
        if x_nodes(j-1) <= x & x <= x_nodes(j) 
            y_eval(i)=(M(j-1)*(x_nodes(j)-x)^3+M(j)*(x-x_nodes(j-1))^3)/(6*(x_nodes(j)-x_nodes(j-1))) ...
                +(y_nodes(j-1)*(x_nodes(j)-x)+y_nodes(j)*(x-x_nodes(j-1)))/(x_nodes(j)-x_nodes(j-1)) ...
                -(x_nodes(j)-x_nodes(j-1))*(M(j-1)*(x_nodes(j)-x)+M(j)*(x-x_nodes(j-1)))/6;
        end
    end
end

function [x, alpha, beta, ier] = tridiag(a,b,c,f,n,iflag)

% function [x, alpha, beta, ier] = tridiag(a,b,c,f,n,iflag)
%
% Solve a tridiagonal linear system M*x=f
%
% INPUT:
% The order of the linear system is given in n.
% The subdiagonal, diagonal, and superdiagonal of M are given
% by the arrays a,b,c, respectively.  More precisely,
%     M(i,i-1) = a(i), i=2,...,n
%     M(i,i)   = b(i), i=1,...,n
%     M(i,i+1) = c(i), i=1,...,n-1
% iflag=0 means that the original matrix M is given as specified above.
% iflag=1 means that the LU factorization of M is already known and is
% stored in a,b,c.  This will have been accomplished by a previous call
% to this routine.  In that case, the vectors alpha and beta should 
% have been substituted for a and b in the calling sequence.
%
% OUTPUT:
% Upon exit, the LU factorization of M is already known and is stored
% in alpha,beta,c.  The solution x is given as well.
% ier=0 means the program was completed satisfactorily.
% ier=1 means that a zero pivot element was encountered and the 
% solution process was abandoned.

a(1) = 0;
if iflag == 0
   % Compute LU factorization of matrix M.
   for j=2:n
      if b(j-1) == 0
         ier = 1;
         return
      end
      a(j) = a(j)/b(j-1);
      b(j) = b(j) - a(j)*c(j-1);
   end
   if b(n) == 0 
      ier = 1;
      return
   end
end

% Compute solution x to M*x = f.
% Do forward substitution to solve lower triangular system.
for j=2:n
   f(j) = f(j) - a(j)*f(j-1);
end

% Do backward substitution to solve upper triangular system.
f(n) = f(n)/b(n);
for j=n-1:-1:1
   f(j) = (f(j) - c(j)*f(j+1))/b(j);
end

% Set output variables.
ier = 0;
x = f;
alpha = a; beta = b;