在Matlab中，我们可以使用pwch()函数来实现
p=pwch(x,y,s)返回插值系数，我们可以使用ppval(p,x0)来计算某一点的插值结果
但是请记住pwch返回的不是每一段的直接多项式的系数，例如下面的表达错误的
>> x = 1:10;
>> y = sin( x );
>> s = cos( x ); % slopes
>> pp = pwch( x, y, s ) ；%interpolation
>>polyval(pp.coefs(1,:),x(1))%这个用法是错误的
而实际上pwch返回的pp.coefs的第i行，是下面多项式的系数
f=coefs(i,1)*(X-breaks(i))^3 +coefs(i,2)*(X-breaks(i))^2 +coefs(i,3)*(X-breaks(i)) + coefs(i,4)
其中，i取从1到length(x)-1
Matlab只提供了怎么计算X对应插值的函数ppval()，没有提供直接多项式系数提取函数。也就是说,你要计算某些点的对应插值时，只需ppval(pp,xx)，这里pp是pwch函数的输出pp形式，必须照写，xx可以是标量或矩阵
而你是想得到形如c(i,1)*(X)^3 +c(i,2)*(X)^2 +c(i,3)*X + co(i,4)的各段上的多项式系数，就要把最上面的多项式 f
下面的代码可以得到你需要的系数：
%by dsj6700417
%see also
%all rights reserved by www.matlabsky.cn
x = 1:10;
y = sin( x );
s = cos( x ); % slopes
pp = pwch( x, y, s ) ;
[breaks,coefs,l,k,d] = unmkpp(pp);
syms X,
c=zeros(9,4);
for i=1length(x)-1)
f=coefs(i,1)*(X-breaks(i))^3 +coefs(i,2)*(X-breaks(i))^2 ...
+coefs(i,3)*(X-breaks(i)) + coefs(i,4);
c(i,:) = sym2poly(f);
end
c
注：如果你要计算x=1时，插值多项式的值，用polyval(pp.coefs(1,:),x(1))是不正确的
按照上面的程序，可以用下面两个办法得到：
>> ppval(pp,1)
ans =
>> polyval(c(1,:),1)
ans =
再注：下面的代码更直接的得到了你想要的系数：
%by dsj6700417
%all rights reserved by www.matlabsky.cn
x = 1:10;
y = sin( x );
z = cos( x );
pp=zeros(length(x)-1,4);
for i=1length(x)-1)
pp %展开后的插值多项式的系数
