% I0.m - modified Bessel function of 1st kind and 0th order.
%
% S = I0(x)
%
% defined only for scalar x >= 0
% based on I0.c

function S = I0(x)

eps = 10^(-9);
n = 1; S = 1; D = 1;

while D > (eps * S),
        T = x / (2*n);
        n = n+1;
        D = D * T^2;
        S = S + D;
end