function p = interpolate(f, xi, x)
	p = zeros(size(x));
	for j = 1:length(x)
		for k = 1:length(xi)
			Lk = 1;
			for l = 1:length(xi)
				if l == k
					continue;
				end
				Lk = Lk * (x(j)-xi(l))/(xi(k)-xi(l));
			end
			p(j) = p(j) + f(xi(k)) * Lk;
		end
	end				
end