clear all
close all

epsilon=12;

[X,Y]=ndgrid(-1:0.01:1);

XY=epsilon*[X(:) Y(:)];
NXY= sqrt(sum(XY.^2,2));
% Z=besselj(0,NXY);  %Bessel s=2
% Z=sqrt(2/pi)*sin(NXY)./NXY; %Bessel s=3
% Z=besselj(1,NXY);  %Bessel s=4
 Z=sqrt(2/pi)*(sin(NXY)-NXY.*cos(NXY))./(NXY.^3); %Bessel s=5

% Z=exp(-NXY.^2); % gaussiana

n=max(size(X));
Zj=reshape(Z,n,n);
H=surf(X,Y,Zj);
%axis square
set(H,'FaceColor','interp','EdgeColor','none')