PROGRAM oscillators
! simulate coupled linear oscillators in one dimension
DIM u(0 to 21),v(0 to 21),usave(0 to 21)
CALL initial(L,u(),v(),t,dt,usave(),mass$,erase$)
DO
   CALL update(L,u(),v(),ke,t,dt,usave())
   CALL animate(L,u(),ke,t,usave(),mass$,erase$)
LOOP until key input
END

SUB initial(L,u(),v(),t,dt,usave(),mass$,erase$)
    DATA 0,0.5,0,0,0,0,0,0,0,0
    DATA 0,0,0,0,0,0,0,0,0,0
    LET t = 0
    LET dt = 0.025
    LET L = 10                    ! number of particles
    SET WINDOW -1,L+1,-4,4
    SET COLOR "red"
    BOX AREA 0.9,1.1,-0.1,0.1
    BOX KEEP 0.9,1.1,-0.1,0.1 in mass$
    SET COLOR "background"
    BOX AREA -0.1,0.1,-0.1,0.1
    BOX KEEP -0.1,0.1,-0.1,0.1 in erase$
    PLOT line -2,0;L+2,0
    FOR j = 1 to L
        READ u(j)                 ! initial displacements
        BOX SHOW mass$ at j-0.1,u(j)-0.1
    NEXT j
    FOR j = 1 to L
        READ v(j)                 ! initial velocities
    NEXT j
    LET u(0) = 0                  ! fixed wall boundary conditions
    LET u(L+1) = 0
    MAT usave = u        ! note use of matrix assignment instruction
END SUB

SUB update(L,u(),v(),ke,t,dt,usave())
    ! Euler-Richardson algorithm
    DIM a(20),amid(20),umid(0 to 21),vmid(20)
    LET ke = 0
    ! K/M equal to unity
    FOR j = 1 to L
        LET usave(j) = u(j)
        LET a(j) = -2*u(j) + u(j+1) + u(j-1)
        LET umid(j) = u(j) + 0.5*v(j)*dt
        LET vmid(j) = v(j) + 0.5*a(j)*dt
    NEXT j
    LET umid(0) = 0
    LET umid(L+1) = 0
    FOR j = 1 to L
        LET amid(j) = -2*umid(j) + umid(j+1) + umid(j-1)
        LET u(j) = u(j) + vmid(j)*dt
        LET v(j) = v(j) + amid(j)*dt
        LET ke = ke + v(j)*v(j)
    NEXT j
    LET t = t + dt
END SUB

SUB animate(L,u(),ke,t,usave(),mass$,erase$)
    LET pe = (u(1) - u(0))^2      ! interaction with left spring
    FOR j = 1 to L
        ! compute potential energy
        LET pe = pe + (u(j+1) - u(j))^2
        ! transverse oscillation
        BOX SHOW erase$ at j-0.1,usave(j)-0.1
        BOX SHOW mass$ at j-0.1,u(j)-0.1
    NEXT j
    PLOT line -2,0;L+2,0
    SET CURSOR 1,1
    SET COLOR "black"
    PRINT using "t = ###.##": t
    LET E = 0.5*(ke + pe)
    PRINT using "E = #.####": E
END SUB