function [fun_val] = neville(x, x_vec, f_vec)

n        = length(x_vec);
Q        = zeros(n);
Q(1:n,1) = f_vec;
for ii = 2:n
    for jj = 2:ii
        Q(ii,jj) = ( ( x - x_vec(ii-jj+1) ) * Q(ii,jj-1) - ( x - x_vec(ii) ) * Q(ii-1,jj-1) );
        Q(ii,jj) = Q(ii,jj) / ( x_vec(ii) - x_vec(ii-jj+1) );
    end
end
fun_val = Q(n,n);