Fourier-Bessel transform

Hi. I would like to calculate tha Fourier-Bessel transform of a known function (a quasi-triangular function). But I am not able to get the result (I expect a sinc^2-like funciton). Could you help to find the errors in the following code? Excuse me, I am a beginner level. Thnaks!!!

 %%%Defocus PSF
  close all;clear all;clc;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
F=8; % f number
f=36; %focal length in mm
L= 600*(10^-6); % mm wavelength
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

r0=linspace(0,5,1501);

T=5;
N=length (r0);
Deltar=r0(2)-r0(1);
Deltaf=(1)/(T);
Fmax=1/(2*Deltar);
nix=linspace(0,Fmax,(Fmax/Deltaf));
u0=L*F*nix;
stepu=u0(2)-u0(1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
MTFZ0(:)=1-(2/pi)*((u0.*((1-u0.^2).^0.5))+asin(u0));

figure (1)
plot (u0,real(MTFZ0),'b','Linewidth',1.5)
grid on;
legend('MTF with Z=0 and Delta=0')
xlabel('optical normalised frequency cycles/mm') 
ylabel('Intensity (Arb. Un.)')

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%PSF calculation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%Z = 0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

taumax0=r0(end)/(L*F);
tau0=linspace(0,taumax0,length(r0));

PSF0=zeros(size(tau0));
for k=1:length(tau0);
    Bes0=zeros(size(u0));
    Integrandus0=zeros(size(u0));
               for i=1:length (u0);
                 Bes0(i)=besselj(0,2*pi*tau0(k)*u0(i));
                 Integrandus0(i)=real(MTFZ0(i)).*Bes0(i);
                 PSF0(k)=2*pi*sum(Integrandus0(i).*u0(i)*stepu);
            end
      end

PSF0=sqrt(PSF0.^2);
rnew0(:)=tau0(:)*(L*F);    

figure (2)
plot (rnew0,PSF0,'b','Linewidth',1.5)
grid on;
legend('PSF0 with Z=0 and Delta=0')
xlabel('radial distance mm') 
ylabel('Intensity (Arb. Un.)')