function [f_approx, coef] = Poly_intp(z_0, Z, f)
%only for z is one dimensional
[~,n] = size(f);
[approx_n,~] = size(z_0);
coef = ones(n,approx_n);
for i = 1:n
    for j = 1:n
        if j ~= i
            coef(i,:) = coef(i,:).*(z_0'-Z(j))/(Z(i)-Z(j));
        end
    end
end
f_approx = f * coef;
