I've run into an obnoxious problem with the below function behavior around zero. Not only does it have problems AT zero due to the divide by zero, it also has problems approaching zero, even for smallish values of nn (eg 11 or 101), which makes patching the divide by zero difficult as I'm not sure how many points surrounding zero to remove.
% fChandra.m  Chandrasekhar function as described in Smith PoP 2005
function fChandra=fChandra(xx) % the Chandrasekhar function used in Dreicer, Phys Rev 1959, eqn 20
fChandra=(erf(xx)xx.*gradient(erf(xx),xx))./(xx.*xx);
fChandra( find(xx(:)==0) )=0.0; % fix divide by zeros
% the function has numerical problems at zero, which the above fix doesn't entirely solve...
%{
nn=1e4+1;
xx=linspace(5,5,nn);find(xx==0)
figure, hold all
plot(xx,fChandra(xx))
plot(xx,gradient(fChandra(xx),xx))
plot(xx,erf(xx))
plot(xx,xx.*gradient(erf(xx),xx))
plot(xx,1./(xx.^2))
legend('chand','grad(chand)','erf','xx*grad(erf,xx)','1/xx^2')
ylim([1.5,1.5])
%}
