Final solution: ``sch.m''

Place the following code block at the end of your ``sch.m'' and run. The code should then start from random functions, perform a finite difference test, orthonormalize, run 400 iterations of your sd(), and compute the final results with your getPsi(). Finally, the code will display the energy of each state along with mesh plots of the values of $\vert\Psi_n(\vec r)\vert^2$ on planes cutting through the center of your cell.

%# 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));
  viewmid(abs(cI(Psi(:,st))).^2,S)
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. Note that s states are spherically symmetric and p states have two lumps separated by a plane of zeros.