import math

def xt(v):
    return (v)

def vt(x, v, c, k, m):
    return (-c*v - k*x)/m

def e_mecanica(x, v, m , k):
    return (m*(v**2)+k*(x**2))/2

m = 2
c = 0.5
k = 1
xi = 0
vi = 0.1
hf = 0.5
h = 0.1
t = 0

k1xt = xt(vi)
k1vt = vt(xi, vi, c, k, m)
k2xt = xt((vi + ((3/4)*k1vt*h)))
k2vt = vt((xi + (3/4)*k1xt*h), (vi + ((3/4)*k1vt*h)), c, k, m)

Em = e_mecanica(xi, vi, m , k) 
t = 0.1
Emnuevo = Em+1 
while (Emnuevo >= (Em/2)):
    xin = xi + ((1/3)*k1xt + (2/3)*k2xt) * h
    vin = vi + ((1/3)*k1vt + (2/3) * k2vt) * h
    k1vt = vt(xin, vin, c, k, m)
    k1xt = xt(vin)
    k2vt = vt((xin + (3/4)*k1xt*h), (vin + ((3/4)*k1vt*h)), c, k, m)
    k2xt = xt( (vin + (3/4)*k1vt*h))
    Emnuevo = e_mecanica(xin, vin, m , k) 
    xi = xin
    vi = vin
    t+=0.1
print('Tiempo: ' + str(t-0.1))