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 = ind1-1;
    end
    ind2 = ind1 + 1;
    yi(i) = y(ind1) + (y(ind2) - y(ind1))*(xi(i)-x(ind1))/(x(ind2)-x(ind1));
end