```python
from sympy import *
```


```python
init_printing()
```


```python
x, dx = var('x, \dot{x}')
m, k, c, g= var('m k c g')
```


```python
T = Matrix([Rational(1, 2) * m * dx**2])
U = Matrix([Rational(1, 2) * k * ((x-1) * Heaviside(x-1))**2 +  - m * g * x])
```


```python
D = T.jacobian(Matrix([dx])).jacobian(Matrix([dx]))
C = T.jacobian(Matrix([dx])).jacobian(Matrix([x]))
G = -(T - U).jacobian(Matrix([x])) 
```


```python
f = Matrix([(- c * dx) * Heaviside(x-1)])
```


```python
p = [m, k, c, g]
v = [x, dx]
```


```python
D_mat = lambdify(p+v, D, modules=['numpy', 'sympy'])
G_mat = lambdify(p+v, G, modules=['numpy', 'sympy'])
f_mat = lambdify(p+v, f, modules=['numpy', 'sympy'])
```


```python
ddq = D.inv().dot(-G+f)
```


```python
ddq_mat = lambdify(v+p, ddq, modules=['numpy', 'sympy'])
```


```python
p1 = [1, 1000, 1000, 9.8]
v1 = [0.0, 0]
```


```python
v2 = [2.0, 0]
```


```python
def foo(t, x, p):
    return np.array([x[1], ddq_mat(*(x+p))])
```


```python
from scipy.integrate import odeint
```


```python
import numpy as np
```


```python
t = np.linspace(0, 10, 100)
```


```python
sol = odeint(foo, t, v1, args=(p1,))
```


```python
y, dy = sol[0, :], sol[1, :]
```


```python
import matplotlib.pyplot as plt
%matplotlib inline
```


```python
plt.plot(y)
```
