function s=cspline(nodes,data,x)

% computes coefficients for natural cubic spline
% at nodes with values data
% evaluate the spline at all points of the array x

N=length(nodes)-1; 
P=length(x); 

% compute the steplength
h=diff(nodes); 
%compute first divided differences of data
D=diff(data)./h;
% compute of D for rhs of linear system
dD3=3*diff(D); 

a=data(1:N)

% generate tridiagonal system 
H=diag(2*(h(2:N)+h(1:N-1))); 
for k=1:N-2
    H(k,k+1)=h(k+1); 
    H(k+1,k)=h(k+1); 
end
c=zeros(1,N+1); 
c(2:N)=H\ dD3'; 

b=D-h.*(c(2:N+1)+2*c(1:N))/3; 
d=(c(2:N+1)-c(1:N))./(3*h); 

% evaluattion 

for i=1:P
    k=1; 
    while (x(i)>nodes(k+1)) & (k<N)
        k=k+1; 
    end
    z=x(i)-nodes(k); 
    s(i)=a(k)+z*(b(k)+ z*(c(k)+z*d(k))); 
end
