def forces(y,t,para):
    nions=y.size/6
    dy=np.empty_like(y)
    ut=np.triu_indices(nions,1)
    wx,wy,b=para
    w=np.array([wx,wy,1.0])
    dy[range(3*nions)]=y[range(3*nions,6*nions)]
    r3=ssd.pdist(y[range(3*nions)].reshape(nions,3,order='F'),'euclidean')**3.0
    for k in range(3):
        X=np.zeros((nions,nions))
        X[ut[0],ut[1]]=(y[ut[0]+k*nions]-y[ut[1]+k*nions])/r3
        X=X-X.T
        dy[range((3+k)*nions,(4+k)*nions)]=np.sum(X,axis=1)-w[k]*y[range(k*nions,(k+1)*nions)]-b*y[range((k+3)*nions,(k+4)*nions)]
    return dy