1;

sample_box=[-2, -2, 2, 2; ...
            inf, 0, 0, inf]; #the classical infinite box

sample_box2=[-2, -2, 2, 2; ...
             0, -3, -3, 0]; # a box in [-2,2] with depth 3

sample_box3=[-2, -2, 2, 2; ...
             2, -3, -3, 2]; # a box in [-2,2] with depth 3

sample_box4=[-1, -1, 1, 1; ...
              0, -6, -6, 0];

sample_box5=[-1, -1, 1, 1; ...
              0,  6, 6, 0]; #this is a wall

function res=set_maxspace(boxes,dx,enlarge=true)
  ##reduce space between points to max dx for plots of solutions
  ##eventually enlarge the interval solution for plots
  xs=boxes(1,1);
  xl=boxes(1,end);
  adj=(xl-xs)*0.1;
  if (enlarge && (isfinite(boxes(2,1))))
    xs -= adj;
    boxes=[[xs;boxes(2,1)],boxes];
  endif
  if (enlarge && (isfinite(boxes(2,end))))
    xl += adj;
    boxes=[boxes,[xl;boxes(2,end)]];
  endif
  n=(xl-xs)/dx+16;
  xx=zeros(1,n);
  yy=zeros(1,n);
  ixx=0;
  ib=1;
  prev_y=boxes(2,1);
  prev_x=boxes(1,1);
  while (ib <= size(boxes,2))
    ixx+=1;
    xb=boxes(1,ib);
    yb=boxes(2,ib);
    if (xb > prev_x + dx*1.1)
      prev_y=prev_y+(yb-prev_y)*dx/(xb-prev_x);
      prev_x=prev_x+dx;
    else
      prev_x=xb;
      prev_y=yb;
      ib+=1;
    endif
    xx(ixx)=prev_x;
    yy(ixx)=prev_y;
  endwhile
  res=[xx(1:ixx);yy(1:ixx)];
endfunction
  
function plot_boxes(boxes)
  ## xs=boxes(1,1);
  ## xl=boxes(1,end);
  ## adj=(xl-xs)*0.1;
  ## xs -= adj; xl += adj;
  ## xx=(0:127)*(xl-xs)/128+xs;
  ## lboxes=[[xs;boxes(2,1)],boxes,[xl;boxes(2,end)]];
  ## yy=interp1(lboxes(1,:),lboxes(2,:),xx);
  ## yy(yy==inf)=NaN;
  ## size(yy)
  ## max(yy)
  ## min(yy)
  xy=set_maxspace(boxes,0.1);
  plot(xy(1,:),xy(2,:));
endfunction

function psi1=schroedinger_eq(v_en,mass,psi0,x)
  psi1=[psi0(2,1),psi0(1,1)*2*mass*v_en(x)];
endfunction

function psij=schroedinger_eq_j(v_en,mass,psi0,x)
  ## the jacobian of the schroedinger eq
  psij=[0,1; ...
        2*mass*v_en(x),0];
endfunction

function [psi,phase_r]=solve_schroedinger1(boxes,energy,mass)
  ##solve a single energy trial schroedinger equation with potential represented by boxes
  n=size(boxes,2);
  psi=zeros(2,n);
  ##psi(1,n) are psi values
  ##psi(2,n) are psi' values
  kk=(boxes(2,1)-energy)*2*mass;
  if (kk < 0)
    ##arbitrary fase, in the continuum, set outgoing wave
    psi(1,1)=1;
    psi(2,1)=-sqrt(kk);
  elseif (kk>sqrt(realmax))
    psi(1,1)=0;
    psi(2,1)=1;
  else
    psi(1,1)=1;
    psi(2,1)=sqrt(kk); ##psi must be asymptotically exp(kx)
  endif
  for j=2:n
    if (boxes(1,j)>boxes(1,j-1))
      x=[boxes(1,j-1),boxes(1,j)];
      psii=[psi(1,j-1),psi(2,j-1)];
      v_en=@(x) segment(boxes(1,j-1),boxes(1,j),boxes(2,j-1),boxes(2,j),x) - energy;
      schr=@(y,x) schroedinger_eq(v_en,mass,y,x);
      schrj=@(y,x) schroedinger_eq_j(v_en,mass,y,x);
      psif=lsode({schr,schrj},psii,x);
      psi(1,j)=psif(2,1);
      psi(2,j)=psif(2,2);
    else
      psi(1,j)=psi(1,j-1);
      psi(2,j)=psi(2,j-1);
    endif
  endfor
  kk=(boxes(2,end)-energy)*2*mass;
  if (kk < 0)
    ##arbitrary fase, in the continuum, set outgoing wave
    psi(1,1)=1;
    psi(2,1)=-sqrt(kk);
    phase_r=psi(1,end)*sqrt(kk)*psi(2,end);
  elseif (kk>sqrt(realmax))
    phase_r=psi(1,end);
  else
    phase_r=sqrt(kk)*psi(1,end)+psi(2,end);
  endif
endfunction

function res=phase_sign(v,en,mass)
  [~,res]=solve_schroedinger1(v,en,mass);
endfunction

function solutions=all_bound_solutions(boxes,mass,top=NaN)
  min_en=min(boxes(2,:));
  if (isnan(top))
    max_en=max(boxes(2,:));
    if (isinf(max_en))
      error("invalid energy interval")
    endif
  else
    max_en=top;
  endif
  max_en=min(max_en,boxes(2,1));
  max_en=min(max_en,boxes(2,end));
  delta=(max_en-min_en)/128;
  en_0=min_en+delta;
  [~,phase_a]=solve_schroedinger1(boxes,en_0,mass);
  v=0;
  for l=2:127
    en_l=min_en+delta*l;
    [~,phase_l]=solve_schroedinger1(boxes,en_l,mass);
    if (phase_a*phase_l<=0)
      ## there must be a solution in this interval, refine it
      en=fzero(@(x) phase_sign(boxes,x,mass),[en_l-delta,en_l]);
      v+=1;
      solutions(v).en=en;
      solutions(v).psi=solve_schroedinger1(boxes,en,mass);
      printf("found a solution at en=%f\n",en);
    else
      printf("no solution in interval [%f,%f]\n",en_l,en_l-delta);
    endif
    phase_a=phase_l;
  endfor
endfunction

function plot_solutions(boxes,psi)
  n=size(psi,2);
  width=size(boxes,2);
  xx=boxes(1,:);
  max_en=realmin;
  min_en=realmax;
  for i=1:n
    psic(i,:)=psi(i).psi(1,:);
    psic(i,:)=psic(i,:)/max(abs(psic(i,:)));
    en=psi(i).en;
    max_en=max(max_en,en);
    min_en=min(min_en,en);
    ## normalize amplitude for plotting
  endfor
  plh=(max_en-min_en)/(n+1); ## try to avoid overlapping picture too much
  plot(xx,boxes(2,:));
  hold on;
  for i=1:n
    en=psi(i).en;
    fmt=sprintf(":;%f;",en);
    plot(xx,psic(i,:)*plh+en,xx,xx*0+en,fmt);
  endfor
  hold off;
endfunction

function box=co_box()
  [r0,e0,k,mass]=CO()
  morse=morse_erk(e0,r0,k);
  x=linspace(1.5,5,20);
  y=morse(x);
  box=[x;y];
endfunction
