function yi = interpolate(x,y,xi)
Ndim = size(y,1);
if any(xi < x(1))  any(xi > x(end))
error('extrapolation needed');
end
yi = zeros(Ndim,length(xi));
for i=1:length(xi)
ind1 = find(x<=xi(i),1,'last');
if ind1 == length(x)
ind1 = ind11;
end
ind2 = ind1 + 1;
yi(:,i) = y(:,ind1) + (y(:,ind2)  y(:,ind1))*(xi(i)x(ind1))/(x(ind2)x(ind1));
end
