import matplotlib.pyplot  as plt
import numpy as np
from scipy.integrate import odeint
Vsat=15.;R1=1.2e3;R2=2.2*R1;alpha=1+R2/R1;Ne=3000;R=1e4;C=1e-7;w0=1/R/C;f0=w0/np.pi/2;Te=25/f0/Ne

def f1(vplus):
    if abs(vplus)<Vsat/alpha:
        return    alpha*vplus
    elif vplus>=Vsat/alpha:
        return    Vsat
    else:
        return    -Vsat

f1vect=np.vectorize(f1)

#u1=np.linspace(-2*Vsat/alpha,2*Vsat/alpha,1000)
#plt.plot(u1,f1vect(u1))
#plt.axis([-2*Vsat/alpha,2*Vsat/alpha,-16,16]);plt.grid()
#plt.show()

def f2(vplus):
    if abs(vplus)<Vsat/alpha:
        return alpha
    else:
        return 0.

def f(y,t_ech):
    return np.array([y[1],(f2(y[0])-3)*w0*y[1]-w0**2*y[0]])

t_ech=np.array([n*Te for n in range(Ne)])
yini=np.array([0.1,0.])
sol=odeint(f,yini,t_ech)
vplus=sol[:,0]
vs=f1vect(vplus)
vmoins=vs/alpha
plt.plot(t_ech,vplus);plt.show()

plt.plot(t_ech,vmoins)
plt.plot(t_ech,vmoins,t_ech,vplus)
plt.plot(t_ech,sol[:,1])
plt.plot(sol[:,0],sol[:,1]);plt.show()