%	Continuous source functions for the diffusion equation
%	Function is normalized with a^2/4pi
clear all
global TAU N
x=[0: 0.02: 5];
tau=[0.2: 0.2: 2];
%	This is tau divided by a^2
integral=tau;
u=x;
for n=1:3
   clf reset
	for i=1:length(tau)
      switch n
      case 1
         u(:)=2*sqrt(tau(i)/pi)*exp(-x(:).^2./(4*tau(i)))...
              	-x(:).*erfc(x(:)./(2*sqrt(tau(i))));
      case 2
         u(:)=(4*pi)^(-1).*expint(x(:).^2./(4*tau(i)));
      case 3
         u(:)=(4*pi.*x(:)).^(-1).*erfc(x(:)./(sqrt(4*tau(i))));
      end
      plot(x,u,'k')
      axis([0 5 0 1])
   	xlabel('x')
   	ylabel('U')
   	title(['Step Response Function for ',num2str(n),'-D Diffusion Equation, tau/a^2 to ',...
      	num2str(tau(end))])
   	hold on
   	TAU=tau(i);
   	N=n;
     	integral(i)=quadl('continf',1.e-6,10);
   end
   N
   integral
   pause
end
