clc; clear; close all;

print_figures = 0;

% Bessel functions of the first kind
m_max = 2000;
dx = 0:0.01:20;
orders = [0, 1, 2];

figure;
for i = 1:length(orders)
    order = orders(i);
    J = BesselFunction(order,m_max);
    plot(dx, polyval(J,dx), 'LineWidth', 2); hold all;
end
axis([0, 20, -0.5, 1.1]);
grid on;
legend('J_0','J_1','J_2','Location','NorthEast');
title('Bessel functions of the first kind for \alpha = 0,1,2');
if print_figures == 1
    filename_string = ['Approximate_Bessel_functions_of_first_type_of_orders_0,1,2'];
    saveas(gcf,filename_string,'pdf');
    saveas(gcf,filename_string,'png');
end
% MATLAB implementation
figure;
for i = 1:length(orders)
    order = orders(i);
    J_values = besselj(order,dx);
    plot(dx, J_values, 'LineWidth', 2); hold all;
end
axis([0, 20, -0.5, 1.1]);
grid on;
legend('J_0','J_1','J_2','Location','NorthEast');
title('Bessel functions of the first kind for \alpha = 0,1,2');
if print_figures == 1
    filename_string = ['MATLAB_Bessel_functions_of_first_type_of_orders_0,1,2'];
    saveas(gcf,filename_string,'pdf');
    saveas(gcf,filename_string,'png');
end


% Bessel functions of the second kind
figure;
for i = 1:length(orders)
    order = orders(i);
    J_values = bessely(order,dx);
    plot(dx, J_values, 'LineWidth', 2); hold all;
end
axis([0, 20, -3.5, 1.5]);
grid on;
legend('Y_0','Y_1','Y_2','Location','NorthEast');
title('Bessel functions of the second kind for \alpha = 0,1,2');
if print_figures == 1
    filename_string = ['MATLAB_Bessel_functions_of_second_type_of_orders_0,1,2'];
    saveas(gcf,filename_string,'pdf');
    saveas(gcf,filename_string,'png');
end