Думаю, если автор ещё не видел вот этой статьи, с ней обязательно стоит ознакомиться. В своё время писал код, опираясь именно на неё.
Разумеется, никаких дифференциальных уравнений.

Плюс, возможно это не сразу становится очевидным: для выполнения действительно точного расчёта порядок действий должен быть сложнее, чем просто выполнить расчёт сразу после обнаружения наложения двух кругов друг на друга, в связи с необходимостью дробить шаг симуляции ∆t.

В статье есть, помимо прочего, расчёт расстояния пробега до столкновения, поэтому, не заостряя внимание на самих формулах, приведу общий план итерации симулятора (для реализации описанного ниже плана потребуется позаимствовать процедуры определения дистанции свободного пробега до столкновения и вычисления векторов скоростей после соударения):

1. просчитать предварительное новое положение каждого из шаров (прибавить ∆t *V_i к текущим координатам каждого i-го шара), и в случае наложения «предсказанных» шариков вычислить для каждого пробег ДО касания (я буду использовать значение доли λ_ij от ∆t (шага симуляции);
   доля шага симуляции (лямбда) — это значение от 0 (почти столкнулись на предыдущем шаге моделирования) до 1 (ещё чуть-чуть, и столкновения на этом шаге бы не было совсем));
    
2. отсортировать пары шаров по полученным значениям долей, чтобы рассматривать сначала те столкновения, которые случились раньше;
    
3. для очередной пары столкнувшихся рассчитать
    
   (а) положения в момент касания (t+∆t*λ_ij),
    
   (б) векторы скоростей (см. статью по ссылке выше),
    
   (в) положения на момент t+∆t (к положениям в момент касания добавить значение вектора скорости после соударения, умноженное на ∆t*(1−λ_ij), корректируя таким образом положение пары на момент t+∆t с учётом столкновения;
    
4. повторять цикл для всех пар шаров, заново определяя наложения после каждого уточнения и начиная расчёт уточнений сначала (для коллизий с самым маленьким значением λ)
    
5. если больше коллизий среди набора предсказанных положений для t+∆t нет, значит, кхм, we could call it a step.

Разумеется, это всё ещё далеко не Virtual pool, но визуально, особенно для моделей с грубым шагом симуляции, это обеспечит ощутимо бóльшую точность. Разумеется, такие заморочки нужны только если таковая нужна. Если нет — из статьи выше достаточно позаимствовать формулы для расчёта углов отскока.

… Помимо этого, при моделировании поведения абсолютно упругих шаров есть тонкий момент, связанный с точным расчётом отскока шара от любого острого угла. Если этот вопрос вдруг также интересен (лично у меня этот вопрос возник непосредственно после решения вопроса с соударениями шаров), могу показать, как считать отскок и для этого случая.