% Example code, % Finite difference solution of Schrodinger equation si_fundconst % fundamental constants alcomp = .3; % aluminum mole fraction d = 10e-9; % quantum well width. c1 = 145e-9; % cladding 1 thickness. c2 = 145e-9; % cladding 2 thickness. d_dep = .3e-6; % depletion width npts = 1024*2; % number of points in mesh % generate mesh; x = linspace(0,c1+c2+d,npts).'; dx = x(2)-x(1); % allocate memory for effective mass and potential meff = zeros(size(x)); V = zeros(size(x)); % Conduction band meffgaas = .067*mo; meffalgaas = (.067+.083*alcomp)*mo; %meffalgaas = (.067)*mo; % using same effective mass % % Heavy hole band % meffgaas = .5*mo; % meffalgaas = (.5+.29*alcomp)*mo; % % % light hole band % meffgaas = .087*mo; % meffalgaas = (.087+.063*alcomp)*mo; % find indicies of cladding, well, cladding indc1 = find(x<c1); indd = find(x>=c1 & x<=c1+d); indc2 = find(x>c1+d); % set effective mass appropriately meff(indc1) = meffalgaas; meff(indd) = meffgaas; meff(indc2) = meffalgaas; % Voltage bias across device if ~exist('Vbias'), Vbias = 3 ;end % Potential due to E-field (eV) Vapplied = linspace(0,Vbias./d_dep*(d+c1+c2),length(V)).'; % Well depth, band discontinuity (eV) % conduction band; V(indd) = -1.247*alcomp*.67; % % valence band; % V(indd) = -1.247*alcomp*.33;% infinite well; % V(indd) = -1e3; % Total potential. V = V + Vapplied; % Set potential to zero at well center, convert to SI units V = (V-mean(V(indd)))*q; % -hbar^2/2/dx^2 common factor dop = ones(size(x))./dx.^2.*(-hbar.^2./2); % plus and minus step effective mass meffp = [meff(2:end);meff(end)]; meffm = [meff(1,:);meff(1:end-1)]; % use either or, to set the effective mass at the boundary %meffp(indc1(end)) = meffp(indc1(end)-1); meffm(indd(1)) = meffm(indd(1)+1); meffp(indd(end)) = meffp(indd(end)-1); %meffm(indc2(1)) = meffm(indc2(1)+1); % upper, lower, and main diagonal, applying flux boundary conditions. udiag = dop(1:end-1)./meffp(1:end-1); cdiag = -dop.*(1./meffp+1./meffm)+V; ldiag = dop(2:end)./meffm(2:end); % generate tri-diagonal matrix A = gallery('tridiag',udiag,cdiag,ldiag); % some options for sparse eigenvalue solver options.issym = 0; options.disp = 0; tic % find 10 eigenvalues of smallest magnitude [PSI,E]=eigs(A,10,'SM',options); toc % These are the energies E = diag(E); % normalize the wavefunction PSI = PSI./(ones(size(PSI,1),1)*sqrt(trapz(x,abs(PSI).^2,1))); % you should plot the PSI's to make sure you have bound states. % figure % doubley2(x,V,x,abs(PSI(:,1:2)).^2); % figure % subplot(211) % plot(x,abs(PSI(:,1)).^2); % subplot(212) %
View Full Document
Unlocking...