%Lagrange 插值
x=[25,40,50,60];
y=[95,75,63,54];
xh=70;% to predict
n = length(x);
m = length(xh);
yh = zeros(1,m);
c1 = ones(n-1,1);
c2 = ones(1,m);
for i=1:n
xp = x([1:i-1 i+1:n]);
yh = yh + y(i)*prod((c1*xh-xp'*c2)./(x(i)-xp'*c2));
end
yh

%牛顿插值
n=10; % number of sample
f=zeros(10,10);
xi=1994:2003;
yi = [67.052,68.008,69.803,72.024,73.400,72.063,74.669,74.487,74.065,76.777];
x=2010;% to predict
for k=1:n
f(k) = yi(k);
end
for i = 2:n
for k = i:n
f(k,i)=(f(k,i-1)-f(k-1,i-1))/(xi(k)-xi(k+1-i));
end
end
disp("差商表");
disp(f);
p=0;
for k = 2:n
t=1;
for j = 1:k-1
t = t*(x-xi(j));
end
p = f(k,k)*t + p;
end
p = f(1,1)+p;
disp("牛顿差值结果为");
p

//newton before
function [P R]=newtoninter(X,Y,x0)
h=abs(X(2)-X(1));
n=find(abs(x0-X)<h);
X=X(n(1):end);
Y=Y(n(1):end);
w=length(X);
R=zeros(w,w);
R(:,1)=Y(:);
for k=2:w;
for j=k:w
R(j,k)=R(j,k-1)-R(j-1,k-1);
end
end
t=(x0-X(1))/h;
T=1;
P=R(1,1);
for m=1:w-1
T=T*(t-m+1);
P=P+R(m+1,m+1)*T/factorial(m);
end
end
%P

//threesimple
xi = -1:0.1:4;
x = 0:3;
y = [3, 5, 4, 1];
pp = csape(x, [1, y, 1], [2, 2]);
yi = ppval(pp, xi);
plot(x, y, 'o', xi, yi, '-');
disp('插值多项式总共分为3段，每个多项式的系数依次为');
for i = 1:size(pp.coefs, 1)
disp(pp.coefs(i, :));
end

https://www.ilovematlab.cn/thread-90910-1-1.html

function f = Hermite(x,y,y_1,x0)
syms t;
f = 0.0;
if(length(x) == length(y))
if(length(y) == length(y_1))
n = length(x);
else
disp('y和y的导数的维数不相等！');
return;
end
else
disp('x和y的维数不相等！');
return;
end
for i=1:n
h = 1.0;
a = 0.0;
for j=1:n
if( j ~= i)
h = h*(t-x(j))^2/((x(i)-x(j))^2);
a = a + 1/(x(i)-x(j));
end
end
f = f + h*((x(i)-t)*(2*a*y(i)-y_1(i))+y(i));
if(i==n)
if(nargin == 4)
f = subs(f,'t',x0);
else
f = vpa(f,6);
end
end
end