function plot_bessel(bchoice)

close all;

if nargin == 0
    bchoice = 'J';
end

% Plot BesselJ
figure(1);
c = {'r','b','m','g','k','y','c'};

switch bchoice
    case 'J'
        x = linspace(0,20,500);
        for n = 0:3
            i = n+1;
            p(i) = plot(x,besselj(n,x),c{i},'linewidth',4-n);
            lstr{i} = sprintf('besselj(%d,x)',n);
            hold on;
        end
        hold on;
        plot(x,0*x,'k');
        ylabel('J0(x)','fontsize',16);
        title('J0(x)','fontsize',18);
        legend(p,lstr,'location','northeast');
        
    case 'I'
        figure(2);
        x = linspace(0,8,500);
        nmaxp1 = 5;
        for n = 0:nmaxp1-1
            i = n+1;
            p(i) = plot(x,besseli(n,x),c{i},'linewidth',nmaxp1-n);
            lstr{i} = sprintf('besseli(%d,x)',n);
            hold on;
        end
        
        axis([0,5,-1,10]);
        ylabel('I0(x)','fontsize',16);
        title('I0(x)','fontsize',18);
        legend(p,lstr,'location','northwest');
        
    case 'E'
        figure(3);
        x = linspace(0,20,500);
        nmaxp1 = 5;
        p(1) = plot(x,besselj(1/2,x),c{1},'linewidth',3);
        lstr{1} = sprintf('besselj(1/2,x)');
        hold on;
        p(2) = plot(x,sqrt(2./(pi*x)).*sin(x),c{2},'linewidth',1);
        lstr{2} = 'sqrt(2/(pi x))*sin(x)';
        hold on;
        
        axis([0,20,-1,1]);
        ylabel('E(x)','fontsize',16);
        title('besselj(n+1/2,x)','fontsize',18);
        legend(p,lstr,'location','northeast');
        
end

set(gca,'fontsize',16);
xlabel('x','fontsize',16);



end