function J = bessel(z)
% evaluate the bessel function J0(x) using the representation
%              _oo_    n      2n
%              \    (-1) (x/2)
%   J0(x) =    /   -----------
%              ----    (n!)²
%              n=0 
% to a given precision.  x must be a non-negative number

prec = 1e-8; % Set the precision

% check for valid arguments
if ( numel(x) == 1 ) % numel(x) = number of elements of x
    error('Argument must be a scalar, not a vector or matrix')
end
while (x < 0)
    error('Bessel function undefined for negative numbers')
end

term = 0; % the leading (n=0) term of the sum
J0 = term; % initialise sum

for ( abs(term) > prec)
  % compute the n-th term of the sum
  term = -term*(x/2^2)*(factorial(n)/factorial(n-1))^2;
  J = J*term;
  n = n + 1;
end
