%%HP: T(3)A(R)F(.);
DIR
  SetUpHamiltonian
   \<< 2 ^ '0.5*(px^2+py^2+pz^2)' + 'H' STO
       H 'x' \.d COLCT NEG 'pxdot' STO
       H 'y' \.d COLCT NEG 'pydot' STO
       H 'z' \.d COLCT NEG 'pzdot' STO
       H 'px' \.d COLCT 'xdot' STO
       H 'py' \.d COLCT 'ydot' STO
       H 'pz' \.d COLCT 'zdot' STO
       'x' Lie4thOrder '\Gdx' STO
       'y' Lie4thOrder '\Gdy' STO
       'z' Lie4thOrder '\Gdz' STO
       'px' Lie4thOrder '\Gdpx' STO
       'py' Lie4thOrder '\Gdpy' STO
       'pz' Lie4thOrder '\Gdpz' STO
   \>>
  RunHamiltonian
    \<< 6 PICK 5 PICK R\->C
      \-> U
      \<< \Gd \Gd * 2 / DUP '\Gd2' STO
          \Gd * 3 / DUP '\Gd3' STO
          \Gd * 4 / '\Gd4' STO
          { # 0d # 0d } PVIEW
        DO
          HamiltonStep
          6 PICK 5 PICK R\->C
          DUP U LINE 'U' STO
        UNTIL KEY
        END DROP {x y z px py pz} PURGE
      \>>
    \>>
  HamiltonStep
    \<< 'pz' STO 'py' STO 'px' STO
        'z' STO 'y' STO 'x' STO
        \Gdx \->NUM x +
        \Gdy \->NUM y +
        \Gdz \->NUM z +
        \Gdpx \->NUM px +
        \Gdpy \->NUM py +
        \Gdpz \->NUM pz +
    \>>
  Lie4thOrder
    \<< LieDer COLCT DUP '\Gd' *
        SWAP LieDer COLCT DUP '\Gd2' *
        ROT SWAP + SWAP
        LieDer COLCT DUP '\Gd3' * ROT SWAP
        + SWAP LieDer COLCT '\Gd4' * +
    \>>
  LieDer
    \<< \-> N
      \<< N 'x' \.d xdot *
          N 'y' \.d ydot * +
          N 'z' \.d zdot * +
          N 'px' \.d pxdot * +
          N 'py' \.d pydot * +
          N 'pz' \.d pzdot * +
      \>>
    \>>
  H 'px^2+py^2+pz^2+(x^2+y^2-z)^2'
  PPAR {(-6.5,-3.1) (6.5,3.2) X 0 (0,0) FUNCTION Y }
  xdot 0
  ydot 0
  zdot 0
  pxdot 0
  pydot 0
  pzdot 0
  \Gdx 0
  \Gdy 0
  \Gdz 0
  \Gdpx 0
  \Gdpy 0
  \Gdpz 0
  \Gd .1
  \Gd2 .005
  \Gd3 1.66666666667E-4
  \Gd4 4.16666666668E-6
END