function yi = interpolate(x,y,xi)
if any(xi < x(1))  any(xi > x(end))
error('extrapolation needed');
end
yi = zeros(size(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
