function f=sourcef(x)
%Integrand for integrating step function over volume
global TAU N
u=x;
fac=x;
switch N
      case 1
         u(:)=2*sqrt(TAU/pi).*exp(-x.^2./(4*TAU))...
              	-x.*erfc(x./(2*sqrt(TAU)));
      case 2
         u(:)=(4*pi)^(-1).*expint(x.^2./(4*TAU));
      case 3
         u(:)=(4*pi.*x).^(-1).*erfc(x./sqrt(4*TAU));
end
switch N
case 1
   fac=1;
case 2
   fac=2*pi.*x;
case 3
   fac=4*pi.*x.^2;
end
f=u.*fac;

