% atom
% One-dimensional confining potential,
% solved through matrix diagonalization
% hbar = m = 1
   
   close;close;close;close;

% Number of grid points: N+1
N = 2000;   % Number of steps between Xmin and Xmax
% if potential number 1 or 2 or 3:
% Xmax = 10; Xmin = -Xmax; %% (symmetric range)

%% if potential number 6:
% Range (from zero to XL, potential infinite for x<0):
XL =100; 
Xmax = XL*(N+1.0)/(N+2.0); Xmin =XL/(N+2.0); 

%  If Potential number 4 or 5 (Box potential without or with sinusoidal bottom), choose 
%  amplitude of the cosine  (v0); and number of periods of the cosine  (nL):
%%%%% v0 = 40; nL=2; Xmin=-nL/2.0*N/(N+2); Xmax=-Xmin;

% Step size:
dx = (Xmax - Xmin)/N;
% Position x is array with values Xmin, Xmin+dx, ... , Xmax:
x = (Xmin):dx:(Xmax);

% Choose type of potential by removing one commmenting % sign:

% potential number 1: Harmonic oscillator:
% V = 0.5*x.^2; 
% potential number 2: Symmetric double well:                     
% V = 0.01*x.^4-0.5*x.^2+6.25;  
% V=x.^4;
% potential number 3: V-shaped potential:      
% V = abs(x);                       
% potential number 4: Ordinary box (use nL=1 above):
% V=0.00*x;
%
% Potential number 5: Box potential with sinusoidal bottom; 
% amplitude of the cosine is v0; number of periods of the cosine is nl:
%
% V = v0*(1+cos(2*pi*(x+nL/2.)));
%
% potential number 6: linear potential, infinite for x<0:
%    V = x;
%
% potential number 7: screened atomic potential, choose l,k and Z; 
% k=100000 gives Coulomb case, no screening, hydrogenlike atom
    l = 0; k = 27; xk = -x./k; Z = 55; xx=x/Z;
V = -(1.0/Z)*(1.0+(Z-1.0).*exp(xk)).*x.^(-1)+l*(l+1).*x.^(-2)./2.0;
fprintf(' N = %5i\n XL = %4.9f\n l = %5i\n k = %5i\n',N,XL,l,k);
% plot screened pot,  Coulomb pot, difference, and effective pots for l=1 and 3
coul = -1.0./x; null=zeros(1,N+1);
% screening correction:
scrcorr = V-coul;
scrcorrplussEb=scrcorr-1/72;
u60=x.*exp(-x/6).*(1-5*x/6+5*(x.^2)/27-5*(x.^3)/(12*27)+x.^4/(24*81)-x.^5/(720*243))/20;
%,(x.*exp(-x/6)*(1-5*x/6+5*(x^2)/27-5*(x^3)/(12*27)+x^4/(24*81)-x^5/(720*243)))^2];
plot(xx,scrcorrplussEb,xx,V,xx,coul,xx,u60,'LineWidth',2);
axis([0 XL/Z -0.1 0.05]);
akse=line(xx,zeros); set(akse,'bk');

break;
% screened effective potentials for l=1 and l=3
Veff1 = V+1./x.^2; Veff3 = V+6./x.^2;
% 
%%figure;
%%plot(x,coul,x,V,x,scrcorr,x,Veff1,x,Veff3,x,null);
%%legend('Coul','screened pot','screening corr.','Veff1','Veff3');
%%xlabel('x=r/a');
%%axis([0 100 -0.3 0.15]);
%%title('potentials and effective potentials');
%%% print -dpng fig42155-11-effpot.png; % creates png file, for linux user
%
% Diagonal elements
d = 1/dx^2+V;
% Off-diagonal elements
e = -1/(2*dx^2);
% Setting up the Hamiltonian matrix, first diagonal terms:
H = diag(d);
% Next, include off-diagonal elements:
H(2:(N+1),1:N) = diag(e*ones(1,N)) + H(2:(N+1),1:N);
H(1:N,2:N+1) = diag(e*ones(1,N)) + H(1:N,2:N+1);
% Diagonalizing the matrix solves the Schrödinger equation.
% The command [S,D] = eig(H) produces matrices of eigenvalues (D) and 
% eigenvectors (S) of matrix H, so that H*S = S*D. Matrix D is the 
% canonical form of H - a diagonal matrix with H's eigenvalues on 
% the main diagonal. Matrix S is the modal matrix - its columns are 
% the eigenvectors of H:
[S,D] = eig(H);
% We store the eigenvalues in the array "eigenvalues":
eigenvalues = diag(D);
% Plot wavefunction for lowest eigenvalue (n = 1):
% n = 1;
% plot(x,(S(:,3)')); 
% Plot absolute square of wavefunction nr 1:
% figure;
% plot(x,(S(:,1)').^2); 
% To display several wavefunctions in one plot:
 figure;
 plot(x,(S(:,1)'),x,(S(:,2)'),x,(S(:,3)'),x,(S(:,4)'),x,(S(:,5)'));
%legend('n=1','n=2','n=3','n=4','n=5');
xlabel('x=r/a');
title('the first five solutions (l=0, no screening, XL=30)');
% print -dpng fig42155-11-3.png   % creates png-file
%%%% plot(x,(1000*S(:,1)'),x,(1000*S(:,2)'),x,(1000*S(:,3)'),x,(1000*S(:,4)'),x,V-3*v0,x,null);
 figure;
plot(x,null,x,(S(:,6)'),x,(5*S(:,6)').^2);
title('solution number six')
% print -dpng fig42155-11-4.png   % creates png-file
%%%%axis([-1 1 -3*v0 90]);
% To display several wavefunctions squared in one plot:
% plot(x,(S(:,1)').^2,x,(S(:,2)').^2,x,(S(:,3)').^2,x,(S(:,4)').^2);
% Plot the potential V(x) in a new figure:
% First make arrays for each of the eigenvalues:
e1=eigenvalues(1)*ones(1,N+1);
e2=eigenvalues(2)*ones(1,N+1);
e3=eigenvalues(3)*ones(1,N+1);
e4=eigenvalues(4)*ones(1,N+1);
e5=eigenvalues(5)*ones(1,N+1);
e6=eigenvalues(6)*ones(1,N+1);
e7=eigenvalues(7)*ones(1,N+1);
e8=eigenvalues(8)*ones(1,N+1);
e9=eigenvalues(9)*ones(1,N+1);
e10=eigenvalues(10)*ones(1,N+1);

% figure;  %plot potential alone
% plot(x,V);
% xlabel('q/L')
% ylabel('V/V_0')
% print -dpng name.png

figure;      % plot of v(x) and the first 10 eigenvalues:
plot(x,V,x,e1,x,e2,x,e3,x,e4,x,e5,x,e6,x,e7,x,e8,x,e9,x,e10);
axis([0 100 -0.3 0.15]);
    % Plot the eigenvalues in a new figure:
    % figure;
% plot(eigenvalues);

% print the first ten eigenvalues, and eigenvsalues multiplied by 2*n², relevant for hydrogenlike atom:
tworydberg=2*13.60569193; hs2overma2=tworydberg*Z^2;
fprintf(' hs2overma2 = %7.8f\n ', hs2overma2);
   for n = 1:10
	     fprintf(' n = %5i   eigenvalue = %4.9f   eigenvalue*2*n² = %4.6f   eigenvalue in eV = %7.6f\n',n,eigenvalues(n),eigenvalues(n)*2*n^2,eigenvalues(n)*hs2overma2)
   end
%eigenvalues(1),eigenvalues(2),eigenvalues(3),eigenvalues(4),eigenvalues(5),eigenvalues(6%),eigenvalues(7),eigenvalues(8),eigenvalues(9),eigenvalues(10)

% print 1/eigenvalues, relevant for hydrogenlike atom:
%  1/eigenvalues(1),1/eigenvalues(2),1/eigenvalues(3),1/eigenvalues(4),1/eigenvalues(5),1/eigenvalues(6),1/eigenvalues(7)
    norm=sum((S(:,6)').^2);
fprintf(' norm = %4.7f\n',norm);
% find radius of psi_n00 orbital (with 90 percent probability inside)
% 
for n=1:5:6
prob = 0.0;
for k=1:N+1
prob = prob + S(k,n).^2;
if abs(prob - 0.9) < 0.001
x90 =  dx*k; fprintf( ' n = %4i   k = %5i  prob = %4.5f   x90 = %4.4f\n',n,k,prob,x90) ;
else end
end % k-loop
end % n-loop
% find radius of unscreened psi_100 orbital (with 90 percent inside), using formula from exercise 6
fr = inline( '(2*r^2+2*r+1)*exp(-2*r)-0.1' ); % pratab p 94
x90_100 = fzero(fr,0.26);                     %   - " -
fprintf(' x90_100 = %3.5f (unscreened)\n',x90_100);
%--------------050406050309000202080109--
