python ball physics simulation

I have seen the great tutorial by Peter Colling Ridge on
http://www.petercollingridge.co.uk/pygame-physics-simulation/
and I am extending the PyParticles script
The code is available on the site(for free), I am using PyParticles4.py

Classes used in the tutorial

The Particle Class
Circular 2d objects with radius,mass,velocity,location
The Spring Class
A spring that binds 2 objects (Particles) and uses the Hooke's law (F = -kx) to determine the interaction between them
The Environment Class
The Environment where the Particles interact

I was wondering if I could to use 2 Particles and make a 'Rod' class (like the Spring class in the tutorial) that had a specific length and didn't allow the particles to come closer go further than that (specified) length.
Also,
Appling a force (when needed) to each Particle such that if one is pulled toward the left, so does the other, but Realistically..
Much like if a 2 different types of balls were joined(from the center) using a steel rod, but in 2-d..
And I don't want to use 3rd party modules

Thanks in advance..

EDIT/UPDATE:
Tried to apply constraint theorem (it failed)
Here's the code:

class Rod:
    def __init__(self, p1, p2, length=50):
        self.p1 = p1
        self.p2 = p2
        self.length = length

    def update(self):
        'Updates The Rod and Particles'
        # Temp store of co-ords of Particles involved
        x1 = self.p1.x
        x2 = self.p2.x
        ###### Same for Y #######
        y1 = self.p1.y
        y2 = self.p2.y

        # Calculation of d1,d2,d3 and final values (x2,y2) 
        # from currently known values(x1,y1)...
        # From Constraint algorithm(see @HristoIliev's comment)
        dx1 = x2 - x1
        dy1 = y2 - y1
        # the d1, d2, d3
        d1 = math.hypot(dx1,dy1)
        d2 = abs(d1)
        d3 = (d2-self.length)/d2
        x1 = x1 + 0.5*d1*d3
        x2 = x2 - 0.5*d1*d3
        y1 = y1 + 0.5*d1*d3
        y2 = y1 - 0.5*d1*d3

        # Reassign next positions
        self.p1.x = x1
        self.p2.x = x2
        ###### Same for Y #######
        self.p1.y = y1
        self.p2.y = y2

A rod in 2D has 3 degrees of freedom (2 velocities/positions + 1 rotation/angular freq).
I would represent the position of the center which is modified by forces in the usual way and calculate the position of the particles using the rotation (for simplicity, about the center of the system) variable.
The rotation is modified by forces by

ang_accel = F * r * sin (angle(F,r)) / (2*M * r^2)

Where

ang_accel is the angular acceleration

F is a force acting on a particular ball so there is 2 torques* that add up as there is two forces that add up (vector-wise) in order to update the position of the center.

r is half of the length
angle(F,r) is the angle between the force vector and the radius vector (from the center to the particle that suffers the force),

So that
F * r * sin (angle(F,r)) is the torque about the center, and
2*M * r^2 is the moment of inertia of the system of two points around the center.