% Archivo Script para Splines cubicos naturales

% Datos
x=[];
y=[];
% Construimos la matriz de coeficientes
n=length(x)
h=diff(x);
ds=h(2:n-2);
dl=h(1:n-2);
dr=h(2:n-1);
dd=2*(dl+dr);
M=diag(dd)+diag(ds,-1)+diag(ds,1);
% Construimos los terminos independientes
q=diff(y)./h;
p=3*diff(q);
% Resolvemos el sistema para obtener los c's
c=p/M;
c=[0 c 0];
% Obtenemos b's y d's
cl=c(1:n-1);
cr=c(2:n);
b=q-h.*(2*cl+cr)/3;
d=(cr-cl)./(3*h);
% Coeficientes del Spline
a=y(1:n-1);
s=[d' cl' b' a'];
% Calculo del Spline
ps=mkpp(x,s);
% Evaluacion y representacion
xx=x(1)::x(n);
yy=ppval(ps,xx);
plot(x,y,'*',xx,yy)