import math

def dx(v):
    return (v)

def dv(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

def RK4(xi,vi,xin,vin,c,k,m,h):
    k1vt = dv(xin, vin, c, k, m)
    k1xt = dx(vin)
    k2vt = dv((xin + (1/2)*h), (vin + ((1/2)*k1vt*h)), c, k, m)
    k2xt = dx( (vin + (1/2)*k1vt*h))
    k3vt = dv((xin + (1/2)*h), (vin + ((1/2)*k2vt*h)), c, k, m)
    k3xt = dx( (vin + (1/2)*k2vt*h))
    k4vt = dv((xin + h), (vin + (k3vt*h)), c, k, m)
    k4xt = dx( (vin + k3vt*h))

    xin = xi + (1/6)*(k1xt + 2*k2xt + 2*k3xt + k4xt) * h
    vin = vi + (1/6)*(k1vt + 2*k2vt + 2*k3vt + k4vt) * h
    return(xin,vin)

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

Em = e_mecanica(xi, vi, m , k) 
Emnuevo = Em
t+=h
while (Emnuevo >= (Em/2)):
    xin,vin=RK4(xi,vi,xin,vin,c,k,m,h)
    Emnuevo = e_mecanica(xin, vin, m , k) 
    xi = xin
    vi = vin
    t+=h
print('Tiempo: ' + str(t-h))