//
// interpolation par krigeage
//
X=[1,2,5,7,8,11,13,16];
Y=[4,6,9,10,8,7,9,5];
deff('y=phi(x)','y=x^3')
n=length(X)
k=2;
C=zeros(n,n)
for i=1:n
    for j=1:n
      C(i,j)=phi(abs(X(i)-X(j)))
    end
end
P=ones(n,k);
for i=1:n
    for j=1:k
      P(i,j)=X(i)^(j-1)
    end
end
A=[C,P;P',zeros(k,k)];
b=zeros(n+k,1);
for i=1:n
  b(i)=Y(i)  
end
z=A\b;
alpha=z(1:n);a=z(n+1:n+k)
t=0:0.01:16;
y=zeros(t)
for i=1:length(t)
    d=0
    for j=1:k
      d=d+a(j)*t(i)^(j-1);
    end    
    w=0;
    for j=1:n
       w=w+alpha(j)*phi(abs(t(i)-X(j))) 
    end
    y(i)=d+w   
end
clf()
plot(t,y)
plot(X,Y,'*')
