用Matlab實現了3次樣條曲線插值的算法。邊界條件取為自然邊界條件，即：兩個端點處的2階導數等于0；

共包含3各個函數文件，主函數所在文件（即使用的時候直接調用的函數）為spline3.m,另外兩個函數文件是在splin3函數文件中被調用的自定義函數。一個是GetParam.m,一個是GetM.m。


%GetParam.m文件的內容：
%根據給定的離散點的橫坐標所構成的向量，計算各個區間段的h值；

function GetParam(Vx,Vy)

global gh;
global gf;
global gu;
global gr;
global gd;
global gff;
global gM;
%global gn;

%n=length(Vx);%length()為向量Vx所含元素的個數；
%n=legth(Vx);
%gn=n;
%n=gn;
n=length(Vx);

gh(1)=Vx(2)-Vx(1);
gf(1)=(Vy(2)-Vy(1))/gh(1);

for i=2:1:n-1%從區間0到區間n-1;
gh(i)=Vx(i+1)-Vx(i);
gf(i)=(Vy(i+1)-Vy(i))/gh(i);

gu(i)=gh(i-1)/(gh(i-1)+gh(i));
gr(i)=1-gu(i);
gff(i)=(gf(i-1)-gf(i))/(Vx(i-1)-Vx(i+1));
gd(i)=6*gff(i);
end

%設置與邊界條件有關的參數；
gM(1)=0;%起點的2階導數；
gM(n)=0;%終點的2階導數；
end


%GetM.m文件的內容：
function GetM(Vx)
global gh;
global gf;
global gu;
global gr;
global gd;
global gff;
global gM;
%global gn;

nn=length(Vx);
%nn=gn;
n=nn-2;

b=zeros(n,1);
A=zeros(n,n);

A(1,1)=2;A(1,2)=gr(2);
b(1)=gd(2)-gu(2)*gM(1);
for i=2:1:n-1
A(i,i)=2;
A(i,i-1)=gu(i+1);
A(i,i+1)=gr(i+1);
b(i)=gd(i+1);
end
A(n,n-1)=gu(n);A(n,n)=2;
b(n)=gd(nn-1)-gr(nn-1)*gM(nn);

X=(inv(A))*b;

for i=2:1:nn-1
gM(i)=X(i-1);
end
end

%主函數文件spline3.m的內容：
function result=spline3(x,Vx,Vy)
global gh;
global gf;
global gu;
global gr;
global gd;
global gff;
global gM;
%global gn;

GetParam(Vx,Vy);
GetM(Vx);

%n=length(Vx);
%n=gn;
n=length(Vx);
nn=length(x);
y=zeros(1,nn);
for j=1:1:nn
i=1;
while(x(j)>Vx(i+1))
i=i+1;
end
sn=i;
t1=(Vx(sn+1)-x(j))^3/(6*gh(sn));
t1=t1*gM(sn);
t2=(x(j)-Vx(sn))^3/(6*gh(sn));
t2=t2*gM(sn+1);
t3=Vy(sn)-gM(i)*((gh(i))^2)/6;
t3=t3*(Vx(sn+1)-x(j))/gh(sn);
t4=Vy(sn+1)-gM(sn+1)*((gh(sn))^2)/6;
t4=t4*(x(j)-Vx(sn))/gh(sn);
y(j)=t1+t2+t3+t4;
end
result=y;
end

函數調用的時候， result=spline3(x,Vx,Vy)，x為代求點的橫坐標向量，
(Vx,Vy)為已知的點的坐標。