function res = eval_poly (x, coeffs, powers)
% Evaluate polynomial at points x (x can be a vector). The
% coefficients of the polynomials are coeffs. The optional paramter powers
% gives the powers that appear in the polynomial. Default is
% [0:length(coeffs)-1].

if size (coeffs, 1) ~= 1 && size (coeffs, 2) ~= 1
    error ('eval_poly:wrongcoeffs', 'Coefficient argument must be a vector.');
end

if size (x, 1) ~= 1 && size (x, 2) ~= 1
    error ('eval_poly:wrongx', 'The argument x must be a vector.');
end

if size (coeffs, 2) ~= 1
    coeffs = coeffs';
end

if size (x, 2) ~= 1
    x = x';
end

n = length(coeffs) - 1; % order of polynomial
if ~exist ('powers', 'var')
    powers = 0:n;
end

res = cell2mat (arrayfun (@(p) (x.^p), powers, 'UniformOutput', false))*coeffs;

end