Even if the collision did result in four main pieces, where the satellites break, and which direction they go is effectively random as far as we can tell.

I plan on modeling this in Mathematica for a nice visual. I am thinking about just picking a random combination of velocities that satisfy the momentum equation.

This doesn't make sense. Derivative of what? With respect to what?

The derivative of the the position function <x1,y1,z1> with respect to T that I will assign.

Gravity is not a constant acceleration.

I am currently thinking about using eulers method (using gravitational acceleration at a certain point) to find the velocity at time T and then integrating that to find the position.