2D Gravity Simualtor (broken)

import math as m
import scipy as SP
import scipy.constants as const
import matplotlib.pyplot as plt
import pylab as PL
from copy import deepcopy as DC
def applyforces(mat,timestep):
    G=const.G
    mat2=DC(mat)
    pi=m.pi
    for i1 in xrange(len(mat2)):
        sv1=0
        sv2=0
        o1=mat[i1]
        x1,y1=o1[0]
        m1=o1[2] #mass in kg
        v1x,v1y=o1[1]
        for i2 in xrange(len(mat2)):
            o2=mat[i2]
            x2,y2=o2[0]
            m2=o2[2] #mass in kg
            v2x,v2y=o2[1]
            d=m.sqrt((x1-x2)**2+(y1-y2)**2) #distance in meters
            if (x1,y1)!=(x2,y2) and m2:
                f=G*((m1*m2)/(d**2)) #grav. force in newton
                a1=f/m1 #acceleration in m/s^2
                a2=f/m2 #acceleration in m/s^2
                v1=a1*timestep #velocity in m/s
                v2=a2*timestep #velocity in m/s
                sv1+=v1
                sv2+=v2
            ang=m.atan2(y1-y2,x1-x2) # angle between points
            v1x,v1y=v1x+sv1*m.sin(ang),v1y+sv1*m.cos(ang)
            v2x,v2y=v2x+sv2*m.sin(ang),v2y+sv2*m.cos(ang)
            dx1=v1x*timestep
            dy1=v1y*timestep
            dx2=v1x*timestep
            dy2=v2x*timestep
            x1,y1=x1-dx1,y1-dy1 #new position
            x2,y2=x2-dx2,y2-dy2 #new position
            mat2[i1][0]=[x1,y1]
            mat2[i2][0]=[x2,y2]
            mat2[i1][1]=[v1x,v1y]
            mat2[i2][1]=[v2x,v2y]
            #elif m2:
            #   mat2[i1][2]+=mat2[i2][2] #add mass
            #   mat2[i1][1][0]+=mat2[i1][1][0] #add vx
            #   mat2[i1][1][1]+=mat2[i1][1][1] #add vy
            #   mat2[i2][2]=0 # remove mass
    return mat2
rx=10
ry=10
t=0
timestep=100 # sec
M=(rx*ry)*10
objmat=[]
scalef=(rx*ry)*5
print "building object matrix..."
for x in xrange(rx):
    for y in xrange(ry):
        mass=M
        vx=0
        vy=0
        objmat.append([[x,y],[vx,vy],mass])
fig=plt.figure()
ax=fig.add_subplot(1,1,1)
t=0
while 1:
    t+=1
    if (t*timestep)%1==0:
        ax.cla()
        for obj in objmat:
            #print "x",obj[0][0],"y",obj[0][1],"vx",obj[1][0],"vy",obj[1][1]
            ax.plot(obj[0][0],obj[0][1],'ro',ms=obj[2]/float(scalef))
        plt.savefig(str(t)+".png")
    print t
    print "applying forces"
    objmat=applyforces(objmat,timestep)
print "done"