PROGRAM sho                       ! simple harmonic oscillator
CALL initial(x,v,w2,t,dt,dt2,nshow)
LET counter = 0
DO
   CALL update(x,v,w2,t,dt,dt2)
   LET counter = counter + 1
   IF counter = nshow then
      CALL output(x,v,t)
      LET counter = 0
   END IF
LOOP until key input
END

SUB initial(x0,v,w2,t,dt,dt2,nshow)
    INPUT prompt "initial position = ": x0  ! meters
    LET t = 0
    LET v = 0                     ! initial velocity
    ! natural (angular) frequency
    INPUT prompt "ratio of k/m = ": w2
    INPUT prompt "time step = ": dt
    LET dt2 = 0.5*dt
    LET show_time = 0.1           ! time interval between output
    LET nshow = int(show_time/dt)
    PRINT tab(7);"time";tab(17);"position";tab(28);"velocity"
    PRINT                         ! skip line
END SUB

SUB update(x,v,w2,t,dt,dt2)
    ! Euler-Richardson algorithm
    LET a = -w2*x
    LET vmid = v + a*dt2          ! velocity at midpoint
    LET xmid = x + v*dt2          ! position at midpoint
    LET amid = -w2*xmid           ! acceleration at midpoint
    LET v = v + amid*dt
    LET x = x + vmid*dt
    LET t = t + dt
END SUB

SUB output(x,v,t)
    ! - print number with leading space or minus sign
    ! % print leading zeroes as '0'
    PRINT using "-----%.####": t,x,v
END SUB