Problem with Solving the Diffusion Equation

[Ross]

I am attempting to simulate the diffusion of implanted particles in a thin film of silicon. I am using the PDE solver to try and do this.

This is the differential equation
function [c,f,s] = Gediffusioneqn(x,t,u,DuDx)
%------------------------- Equation for the Diffusion coefficient
T = 900 ; % Temperature in Celsius
E1Ge = 5.08; % eV
Kb = 8.6173324*(10^-5); %Boltzmann constant
DGe = 7.55*(10^3)*exp(-E1Ge/(Kb*(T+273.15))) ; %diffusion coefficient, cm^2/s
%-------------------------
c=1;
f=DGe*DuDx; %This is the diffusion equation
s=0;

These are the boundary conditions they are supposed to represent that the particles cannot diffuse out of the silicon
function [pl,ql,pr,qr] = bcge(xl,ul,xr,ur,t)
pl = 0;
ql = 1;
pr = 0;
qr = 1;

These are the initial conditions
function value = Geinitial(x) ;
q = 5e16; Dose of implant
T = 900 ; % Temperature
E1Ge = 5.08; % eV
Kb = 8.6173324*(10^-5); %Boltzmann constant
DGe = 7.55*(10^3)*exp(-E1Ge/(Kb*(T+273.15))) ; %diffusion coefficient, cm^2/s
value = (q/(2*(pi*DGe)^(1/2)))*exp(((x-50*10^-7)^2)/(4*DGe)); % This is the %generic solution for an implanted film diffusing through an infinite medium with the % time set to 1 second

Code for diffusion solver
%Mesh
x = linspace (0,220*10^(-7),2201);
t = linspace (0,60,61);
m = 0;
sol = pdepe (m,@Gediffusioneqn,@Geinitial,@bcge,x,t);
u = sol(:,:,1);
plot(x,u(60,:))


When I run this I get the following warning/error message when I run my matlab code
Warning: Failure at t=0.000000e+00. Unable to meet integration tolerances without reducing the step size below the
smallest value allowed (7.905050e-323) at time t.
> In ode15s at 669
In pdepe at 317
In Gediffusionsolver at 5
Warning: Time integration has failed. Solution is available at requested time points up to t=0.000000e+00.
> In pdepe at 323
In Gediffusionsolver at 5
Index exceeds matrix dimensions.

Error in Gediffusionsolver (line 7)
plot(x,u(60,:)) .
]

Any help would be greatly appreciated.

Cheers,
Ross

[Bill Greene]

"Ross" wrote in message <likmbc$3ke$1...@newscl01ah.mathworks.com>...

The first problem I notice is that your initial conditions are infinite at all points.
Just do:

Geinitial(x(1))

There is no reason to call pdepe until you get this problem resolved.

Beyond that, here are a couple more observations:

You don't need to do something like 50*10^-7). It is much simpler to
write this as 50e-7. The values are the same; its just that the second is
easier to read.

You have 2201 mesh points. This is an extremely fine mesh. I doubt if
such a fine mesh is required but, even if it is, I suggest starting with something
much coarser (say 20) until you get all other problems resolved. Then you can
refine the mesh somewhat until the solution is sufficiently accurate.

Bill