Earlier versions of my physics engine used a penalty-based method for collisions. While simple to implement, it resulted in "spongy" interactions. To fix this, I refactored the solver to use Impulse Resolution.
The Math
The magnitude of the impulse scalar $j$ is calculated using the coefficient of restitution $e$ (bounciness) and the relative velocity along the normal.
j =
-(1 + e) ċ (vrel ċ n)
mA-1 + mB-1
Implementation
Here is the simplified C# implementation used in the solver loop.