Не очень понятен уровень детализации.
1. Задача гравитационного взаимодействия программируется сравнительно просто
ru.wikipedia.org/wiki/Задача_трёх_тел
у вас будет не три тела, а много (N тел). Т.е. вам нужно будет решать систему из 3N обыкновенных дифуравнений.
Решение будет неправильным (по причине неустойчивости), но в общем похожим на правду.
2. Также вам надо учесть, что у вас не материальные точки, а шары - т.е. при движении звезд учитывать расстояние не только между центрами, но и между поверхностями - когда поверхности соприкоснулись - начинается сливание.
3. Как описать сливание звезд. Тут уже сложно. Звезда - это жидкость. Если хотите видеть как меняется форма поверхности - нужно решать уравнение Навье-Стокса (ну или возможно просто Эйлера) с введенными электромагнитными и гравитационными силами и отслеживать форму поверхности каким-либо методом - например VOF (Volume of fluid). ну это так как на самом деле. Если просто поиграться - то зная суммарный объем можно просто слить звезды как капельки не думая о физике.
Успехов!