anchor.diff:    close(unit=10)
   
     !---------------------------------------------
  -  ! Compute initial speeds
  +  ! Compute initial speeds, accelerations
     !---------------------------------------------
     call initial_speeds(r, v)
  +  call compute_force(m, r, a)
   
     !---------------------------------------------
     ! Loop over time
  
positive.code: 
    real(kind = xp), dimension(:), allocatable :: m
    real(kind = xp), dimension(:, :), allocatable :: r, v, a
    integer :: iter = 0
    integer :: dump_freq = 1
    integer :: maxiter = 10000
  
negative.code:      !-------------------------
       ! Compute acceleration
       !-------------------------
       call integrate(v, a, dt/2)
       call integrate(r, v, dt)
       call compute_force(m, r, a)
       call integrate(v, a, dt/2)
       call compute_energy(m, r, v)
    end do
  
       if (mod(dump_freq, iter) == 0) then
          call write_dump(5, r, v, npoints)
  