%1.5 harmoniccode.m
%The parameters are L=10, n=200 and k=20, which means that the
%interval is [-10,10], there are 200 meshpoints in [0,1] and we want
%the first 20 eigenvalues.
% to plot the wave function of the jth eigen state:
% plot(x,V(1:dim,k-j+1));
% We can convince ourselves that [-10,10] is big enough by looking at the
% graph of the highest wave function. SEE harmonic.mat****
function harmonic(L,n,k) % positive integers L,n,k
h=1/n; N=2*n*L+1; % N points in interval [-L,L]
K=toeplitz([2 -1 zeros(1,N-2]); % second difference matrix
V=diag((-L:h:L).^2/2); % diagonal matrix from x^2/2
H=K/(2*h^2)+V; % Hamiltonian = kinetic + potential
[V,F]=eig(H,k,'sm'); % k eigenvectors and eigenvalues
E=diag(F) % smallest k eigenvalues of H
j=1:k;
plot(k+1-j,E); % MAYBE plot(1:k,E)?