>>716 にあげた式を元に数値計算
だいたい減衰振動的な挙動になりました。
振動周期と緩和時間は何か簡単な関係式で表せるのかもしれません。

(たぶんインデント崩れるだろうから実行時は適宜補う事)
import numpy as np
from scipy.integrate.odepack import odeint
import pylab as plt

# x = [x, v, E1, E2], V1+V2=2, E1+E2=2, N1=N2=1
def func(x, t, m, k):
return[ x[1],
(x[2]/x[0]-x[3]/(2-x[0]))/m,
+k*(x[3]-x[2]) -x[2]/x[0]*x[1],
-k*(x[3]-x[2]) +x[3]/(2-x[0])*x[1] ]

t = np.arange(0, 40, 0.01)
m = 1.5
k = 2.0
x_ini = [0.5, 0, 0.2, 1.8]
x = odeint(func, x_ini, t, args=(m,k))

# V1
plt.plot( x[:,0] )
plt.show()

# E1
plt.plot( x[:,2] )
plt.show()

# E3
plt.plot( 2-x[:,2]-x[:,3] )
plt.show()