Hello MB,
[1] Dimensionally, k can't be wavelength. Maybe you meant the wave number, 2pi / lambda.
[2] The integrand increases proportionally to r^2, which is certainly possible, but somewhat suspicious. Are you certain that an r^2 factor for the volume element has not already been inserted into ( kr j_n(kr) )^2 as the factors of r there? That happens sometimes.
You can use the identity
d/dz Jn(z) = -J(n+1,z) + (n/z)*Jn(z)
along with the identity you have already
jn(z) = sqrt((pi/2)/z) besselj(n+1/2,z)
to obtain
d/d[kr] (kr j_n(kr)) = -kr j_n+1(kr) + (n+1) j_n(kr)
The following function is the integrand, which can be plugged into the 'integral' function, integrating from a to b:
function y = fun(r)
k = 2.2; n = 3; % for example
kr = k*r;
y = ((-kr.*js(n+1,kr) + (n+1)*js(n,kr))).^2.*r.^2;
function y = js(n,x) % nested function
% spherical bessel function j, first kind
% y = js(n,x);
y = sqrt((pi/2)./x).*besselj(n+1/2,x);
% get rid of nans
if n==0
y(x==0) = 1;
elseif n >0
y(x==0) = 0;
else
y(x==0) = (-1)^(n-1)*inf;
end
end % end js
end % end fun
