function y = mybessel(x)
%
%   Computes the 0th-order Bessel function of the first kind
%

K = 32;

bessel_coef =  zeros(1,K);

kfac = 1;
toggle = 1;
two_to_k = 1;

for k = 1:K
    toggle = -toggle;
    kfac = kfac*k;
    two_to_k = two_to_k * 2;
    bessel_coef(k) = toggle/(kfac*two_to_k)^2;      % compute power series coefficients in advance
end

x = x.^2;

    y = x .* bessel_coef(K);
for k = K-1:-1:1
    y = x .* (bessel_coef(k) + y);          % Horner's method
end

y = 1 + y;

end

x = linspace (-14, 14, 28*4096+1);

J_0 = besselj(0, x);            % MATLAB's bessel

I_0 = mybessel(x);              % my bessel

figure(1)
plot(x, J_0)
hold on
plot(x, I_0)
hold off

figure(2)
plot(x, I_0 - J_0)              % plot error