PROGRAM eigen
CALL parameters(V0,a,xmax,dx)
CALL plot_potential(V0,a,xmax)
CALL Euler(V0,a,dx,xmax)
END

SUB parameters(V0,a,xmax,dx)
    INPUT prompt "magnitude of well depth = ": V0
    INPUT prompt "half width of well = ": a
    INPUT prompt "step size = ": dx
    INPUT prompt "maximum value of x to be plotted = ": xmax
END SUB

SUB plot_potential(V0,a,xmax)
    ! draw potential well
    OPEN #1: screen 0,1,0,0.4
    SET COLOR "red"
    SET WINDOW -1.05*xmax,1.05*xmax,-1.1*V0,1.1*V0
    PLOT LINES: -xmax,0;-a,0;     ! horizontal line
    PLOT LINES: -a,0;-a,-V0       ! vertical line
    PLOT LINES: -a,-V0;a,-V0;a,-V0;a,0;a,0;xmax,0
END SUB

SUB Euler(V0,a,dx,xmax)
    DECLARE DEF V
    OPEN #2: screen 0,1,0.4,1
    SET WINDOW -1.01*xmax,1.01*xmax,-4,4
    INPUT prompt "even or odd parity (1 or -1) = ": parity
    CLEAR
    DO
       INPUT prompt "E = ": E
       IF E = 0 then EXIT SUB
       SET COLOR "blue"
       IF parity = -1 then
          LET phi = 0             ! initial values at x = 0
          LET dphi = 1            ! first derivative
       ELSE
          LET phi = 1
          LET dphi = 0
       END IF
       LET x = 0
       DO                         ! compute wave function
          LET x_old = x
          LET phi_old = phi
          LET x = x + dx
          LET d2phi = 2*(V(x,V0,a) - E)*phi      ! dimensionless units
          LET dphi = dphi + d2phi*dx   ! Euler-Cromer algorithm
          LET phi = phi + dphi*dx
          ! plot wave function
          PLOT LINES: x_old,phi_old; x,phi
          PLOT LINES: -x_old,phi_old*parity; -x,phi*parity
       LOOP until x > xmax
       SET COLOR "black"
    LOOP
END SUB

DEF V(x,V0,a)                     ! potential function V(x)
    IF abs(x) > a then
       LET V = V0
    ELSE
       LET V = 0
    END IF
END DEF