% 分段三次埃尔米特插值
% 注释的地方需要修改
clear
xd=[1,2,1]; % 插值节点，逗号分隔
yd=[2,3,3]; % 插值节点对应的函数值，逗号分隔
ydot=[0,-1,3] % 插值节点对应的导数值，逗号分隔
x0=1.7; % 求值节点
n=length(xd);
h=diff(xd);
pp=zeros(n-1,6);
y=@(x)zeros(size(x));
for k=1:n-1
pp(k,1:4)=2*yd(k)/h(k)^3*poly([xd(k+1),xd(k+1),(3*xd(k)-xd(k+1))/2])+...
-2*yd(k+1)/h(k)^3*poly([xd(k),xd(k),(3*xd(k+1)-xd(k))/2])+...
ydot(k)/h(k)^2*poly([xd(k),xd(k+1),xd(k+1)])+...
ydot(k+1)/h(k)^2*poly([xd(k),xd(k),xd(k+1)]);
pp(k,5:6)=xd(k:k+1);
y=@(x)y(x)+polyval(pp(k,1:4),x).*(x>=xd(k) & x<xd(k+1));
end
y0=feval(y,x0);
y0(x0==xd(end))=yd(end);
vpa(y0,7) % 结果有效数字位数