用matlab写一个三次样条插值的代码

function [S] = cubicSpline(x,y)
%三次样条插值
%输入x,y分别是n+1个节点和对应的函数值，其中x必须是递增的
%函数返回一个三次样条,S(x(k))=y(k),k=1:n+1
%
%Author:Jianhua Huang(285183443@qq.com)
%Date:2019/5/27

n=length(x)-1;

%拟合三次函数,得到矩阵A,b
A=zeros(4*n,4*n);
b=zeros(4*n,1);
row=1; %行指针
for i=1:n
%Si(x(i))=y(i),Si(x(i+1))=y(i+1),Si'''(x(i))=Si'''(x(i+1))=0
A(row,4*(i-1)+1:4*i)=[x(i)^3 x(i)^2 x(i) 1];
b(row)=y(i);
row=row+1;
A(row,4*(i-1)+1:4*i)=[x(i+1)^3 x(i+1)^2 x(i+1) 1];
b(row)=y(i+1);
row=row+1;
if i~=n
%Si-1(x(i+1))=y(i+1);
A(row,4*(i-1)+1:4*i)=[x(i+1)^3 x(i+1)^2 x(i+1) 1];
A(row,4*i+1:4*(i+1))=[-x(i+1)^3 -x(i+1)^2 -x(i+1) -1];
A(row,2*i+1:2*i+2)=[2*x(i+1) 1];
row=row+1;
%Si-1''(x(i+1))=Si''(x(i+1));
A(row,4*(i-1)+1:4*i)=[6*x(i+1)^2 2*x(i+1) 0 0];
A(row,4*i+1:4*(i+1))=[-6*x(i+1)^2 -2*x(i+1) 0 0];
A(row,2*i+1:2*i+2)=[2 0];
row=row+1;
end
end

%解线性方程组
c=A\b;

%构造三次样条
S=cell(1,n);
for i=1:n
S{i}.a=c(4*(i-1)+1);
S{i}.b=c(4*(i-1)+2);
S{i}.c=c(4*(i-1)+3);
S{i}.d=c(4*(i-1)+4);
S{i}.x=x(i:i+1);
end

end

%插值函数
function [y] = cubicSplineInterp(S,x)
%三次样条插值函数

n=length(S);
y=zeros(size(x));
for i=1:length(x)
for j=1:n
if x(i)>=S{j}.x(1) && x(i)<=S{j}.x(2)
y(i)=S{j}.a*(x(i)-S{j}.x(1))^3+...
S{j}.b*(x(i)-S{j}.x(1))^2+...
S{j}.c*(x(i)-S{j}.x(1))+...
S{j}.d;
break;
end
end
end

end