
%  The script file table3.m computes and plots the real spectrum
%  of the QES potential p(z) = -z^4-2bz^2-6iz.   It reproduces 
%  Figure 3 of the paper:
%  JAC Weideman, "Spectral Differentiation Matrices for the Numerical
%  Solution of Schroedinger's equation" (submitted).

hold off;
N = 65; mu = 1; c = 5;  alpha = pi/6;  % Good parameters determined empirically
 
[w, DM] = herdif(N, 2, c);   % Assemble Hermite DM of size (N+1)x(N+1) 
DH1 = DM(:,:,1);             % and extract first
DH2 = DM(:,:,2);             % and second derivatives.  

     z = i*mu*(sin(alpha)+sin(i*w-alpha));   % Define hyperbola
  dwdz = -1/mu*sec(i*w-alpha);               % and change variables.
d2wdz2 = i/mu^2*sec(i*w-alpha).^2.*tan(i*w-alpha);

D2 = diag((dwdz).^2)*DH2+diag(d2wdz2)*DH1;   % Assemble second derivative
                                             % in new coordinate system              

for b = linspace(-4,4,1000);    % Start loop over b. 
 
    A = -D2+diag(-z.^4-2*b*z.^2-6*i*z);
    E = eig(A);

ell1 = find(abs(imag(E)) < 1e-5);  % Retain only real and small e-values 
  Ereal= E(ell1);
ell2 = find(abs(Ereal) < 31);
  ERsmall = Ereal(ell2);
  
  plot(b*ones(size(ERsmall)),ERsmall,'.','MarkerSize',12); hold on;   % Plot
end

axis([-4 4 -5 30])
xlabel('b','FontSize',16)
ylabel('E','FontSize',16)
plot([0 0],[-5 30]); 
plot([-4 4],[0 0]);



