Final solution: ``sch.m''

Place the following code block at the end of your ``sch.m'' and run. (Note: if you wish to regain control before all of the iterations are complete, you have the option of hitting Ctrl-C once (and then $<$enter$>$ if there is no immediate response. Be patient: if you hit Ctrl-C more than once, octave will likely crash, destroying your octave session and creating a large octave-core file.)

The code block below starts from random functions, performs a finite difference test, orthonormalizes the starting functions, runs 400 iterations of your sd(), and computes the final results with your getPsi(). Finally, the code displays the energy of each state along with grey-scale density plots of the values of $\vert\Psi_n(\vec r)\vert^2$ on planes cutting through the center of your cell. For the graphics to function, be sure to copy ``smooth.m'' and ``ppm.m'' from Appendices C and D into your octave directory.

%# Allow for more digits in printouts
format long

%# Converge using steepest descents
W=sd(W,400);

%# Extract and display final results
[Psi, epsilon]=getPsi(W);

for st=1:Ns
  printf('=== State # %d, Energy = %f ===\n',st,epsilon(st));
  dat=abs(cI(Psi(:,st))).^2;
  for k=1:3
    sl=slice(dat,S,S(k)/2,k);
    name=sprintf("psi%d_m%d.ppm",st,k);
    ppm(name,sl*0.3,sl,sl); system(["display " name "&"]);
  end
end

For your reference, with an angular frequency of the oscillator of $\omega=2$, the analytic result is that the first lowest four states have energies 3, 5, 5 and 5, with spatial character s, p, p and p, respectively. To view all of the electron density slices, place the mouse over the graphics window. (You may have to click on the title bar to get the graphics window's ``attention''.) Then, to flip through the pictures, hit the space bar while the mouse is over the image.