Spline Natural en MatLab

function [H,Y,h]=spline_natural(x,y)
% Syntaxis: [yy]=spline_cubico(x,y,xx,tipo,pend)
% ----------------------------------------------------------------------
% Devuelve el valor interpolado yyen un determinado punto xx solicitado
% mediante el empleo de un Spline Cubico de tipo "natural", "anclado" o
% "periodico".
%
% Entrada:
% x --> Vector que contiene las abcisas de losnodos de interpolacion
% y --> Vector que contiene las ordenadas de los nodos de interpolacion
% xx --> Abcisa del nodo que se desea incorporar
% tipo --> Tipo de Spline que se desea emplear% ["natural", "anclado", "periodico" {default | "natural"}]
% pend --> Vector que contiene las pendientes en los extremos para el caso
% del Spline "anclado"
%
%
%Salida:
% yy --> Valor interpolado para el nodo con abcisa xx
%
%
% Ejemplo:
%
% x=[1; 2; 4; 5; 8]; y=[2; 3; 1; 4; 9]; xx=2.5;
% [yy]=spline_cubico(x,y,xx,'natural')
%
% yy =
% 2.4735%
%
% Ver tambien: interpola_lagrange, interpola_newton
%
% ----------------------------------------------------------------------
%
% Elaborado por:
%
% Msc. Alexeis Comanioni Guerra
%Lic. Vianka Orovio Cobo
%
% Revision: 1.0
% Fecha: 30/09/2007


%% Determinacion de Spline Cubico Natural


% Conformacion de la matriz H
nx = length(x);
h = diff(x);

d1 =h(2:end-1)/6;
d2 = (1/3)*( h(1:end-1)+h(2:end) );

H = full( gallery('tridiag',d1,d2,d1) );
nH=nx-2;
% Conformacion del vector columna Y
Yt = diff(y)./h;
Y = Yt(2:end)-Yt(1:end-1);

% Resoluciondel SEL (metodo de gauss) y obtencion del vector M
%M = zeros(nH+2,1);
%M(2:nH+1,1)=gauss(H,Y);

% Determinando el valor interpolado
%i = nx;
%while x(i) > xx
% i = i-1;
%end

%u =xx-x(i); v = x(i+1)-xx;
%yy1 = (1/6)*(1/h(i))*( M(i)*v^3 + M(i+1)*u^3 );
%yy2 = (1/h(i))*( v*y(i)+u*y(i+1) );
%yy3 = (h(i)/6)*( M(i)*v + M(i+1)*u );

%yy = yy1 + yy2 - yy3;