% Fourier Bessel series 2
clear				% Clear previous definitions
figure(101)         % Assign figure number
clf                 % Clear previous figures
hold off            % Start with fresh graph

NptsX=1000;         % number of x pts             
x = linspace(0,4,NptsX);
g=@(z) besselj(0,z);
bz(1) = fsolve(g,2.4);
lam(1)=bz(1)/4;     % sqrt(e.v.)
evlam(1)=lam(1)^2;  % 1st e.v.
for n=2:50
    bz(n)=fsolve(g,bz(n-1)+pi);
    lam(n)=bz(n)/4;     % sqrt(e.v.s)
    evlam(n)=lam(n)^2;  % n^{th} e.v.
end
f1 = zeros(1,NptsX);
for n=1:50
    lamn = lam(n);
    fdenom = @(x) x.*(besselj(0,lamn*x)).^2;
    bdenom(n) = integral(fdenom,0,4);
    fnum = @(x) 4*x.*(besselj(0,lamn*x));
    bnum(n) = integral(fnum,1,3);
    b(n) = bnum(n)/bdenom(n);    % Fourier coefficients
    fn = b(n)*besselj(0,lamn*x); % Fourier function(n)
    f1 = f1+fn;
end
plot(x,f1,'b-','LineWidth',1.5);
hold on
plot([1 3],[4 4],'k-','LineWidth',1.5);
plot([0 1],[0 0],'k-','LineWidth',1.5);
plot([3 4],[0 0],'k-','LineWidth',1.5);
plot([1 3],[2 2],'ko','MarkerSize',5,'MarkerFaceColor','k');
grid;
h = legend('Fourier Bessel series, $n = 50$',...
     '$f(x)$', 'Location','south');
set(h,'Interpreter','latex')
h.FontSize = 10;
xlim([0,4]);
ylim([-1 5]);
xlabel('$x$','FontSize',12,'interpreter','latex');  
ylabel('$y$','FontSize',12,'interpreter','latex');
set(gca,'FontSize',12);         % Axis tick font size
print -depsc bessel_fourier2.eps
