% RBF of high order fundamental solution of 3D Helmholtz operator

function hhelm3d=hhelm3d(r,Lamda,nj)

  A=-Lamda/4/pi;

 if(nj==0&r==0)
     hhelm3d=A;
 elseif(nj==0&r~=0)
     hhelm3d=A*sqrt(pi/2)*(Lamda*r)^(-1/2)*BesselJ(1/2,Lamda*r);
 else
     for j=1:nj
         A=-A/(2*j*Lamda^2);
     end
         hhelm3d=A*sqrt(pi/2)*(Lamda*r)^(nj-1/2)*BesselJ(nj+1/2,Lamda*r);
 end