%program for computing the internal potential in a semiconductor

l = 10; n = 80;
u0 = 0.2                      %external potential
xl = (0:l/n:l)'               %mesh
cl = tanh(20.0.*(xl./l-0.5))  %concentration of impurities
phid = asinh(cl*0.5)          %builtin potential
h = xl(2:n+1)-xl(1:n)         %steps
b = 1 ./h
                              %tridiagonal matrix A 
A = spdiags([b(2:n) -(b(1:n-1)+b(2:n)) b(1:n-1)],-1:1,n-1,n-1)

hs = (h(1:n-1)+h(2:n))*0.5    %average steps
c = cl(2:n)

f0 = -hs.*c
f0(1) = f0(1)-(phid(1)+u0)*b(1); 
f0(n-1) = f0(n-1)-phid(n+1)*b(n); 
phi =  phid(2:n)              %initial approximation

nor = 1                       %iteration
while nor > 0.001 
   ex = exp(phi); f =  f0+hs.*(ex-1./ex)-A*phi
   D = spdiags(hs.*(ex+1./ex),0,n-1,n-1)
   H = (A-D)\f
   nor = norm(H,inf)
   phi = phi+H
end
phiv = [phid(1)+u0; phi; phid(n+1)] %with boundary conditions

clf;
plot(xl,phiv,'-r', xl,phid,':b')
text(4.6,-0.4,'Builtin potential');
text(6,0.3,'Internal potential');
xlabel('x'); ylabel('phi');        
                              %end

