Identification Division.
       Program-Id. Dipole.
      * from Stauffer/Stanley
      * From Newton to Mandelbrot
      * Springer 1996
       Data Division.
       Working-Storage Section.
         01  occurs 20.
           05 ph        pic S9(4)V9(4) occurs 20.
         77 vx          pic S9(4)V9(4).
         77 vy          pic S9(4)V9(4).
         77 dt          pic S9(4)V9(4).
       Procedure Division.
        call hgr hcolor:7
        call hplot using 130,1,to,130,159
        call hplot using 1,80,to,260,80
        L10.
        display "x,y="
        accept x,y
        move 0.1 to dr
        move 10 to sc
        L30.
        compute r2=x*x+y*y
        compute r5=function pow(r2,2.5)
        compute ex=(3*x*x-r2)/r5
        compute ey=3*y*x/r5
        compute x=x+ex*dr*r2
        compute y=y+ey*dr*r2
        display x, y
        call hplot using 130+sc*x,80+sc*y
        call hplot using 130-sc*x,80+sc*y
        call hplot using 130+sc*x,80-sc*y
        call hplot using 130-sc*x,80-sc*y
        if x>0 then go to L30
        goto L10.
       End-Program Dipole.