#4.1.7
a=1
t=srange (0, 50, 0.1)
var ("X", "V")
sol= desolve_odeint((V, -X-(a*V^3-V)), ics= [5, 10], dvars=[X, V], times=t)
list_plot(list(zip(t, sol[:,0])), axes_labels=["time", "Vprime"]) + list_plot(list(zip(t, sol[:,1])), color="red")

#4.1.7
a=1
t=srange (0, 50, 0.1)
var ("X", "V")
sol= desolve_odeint((V, -X-(a*V^3-V)), ics= [5, 10], dvars=[X, V], times=t)
list_plot(list(zip(sol[:,0], sol[:,1])), axes_labels=["Xprime", "Vprime"], plotjoined=True)

#4.1.7
a=5
t=srange (0, 50, 0.1)
var ("X", "V")
sol= desolve_odeint((V, -X-(a*V^3-V)), ics= [5, 10], dvars=[X, V], times=t)
list_plot(list(zip(t, sol[:,0])), axes_labels=["time", "Vprime"]) + list_plot(list(zip(t, sol[:,1])), color="red")

#4.1.7
a=5
t=srange (0, 50, 0.1)
var ("X", "V")
sol= desolve_odeint((V, -X-(a*V^3-V)), ics= [5, 10], dvars=[X, V], times=t)
list_plot(list(zip(sol[:,0], sol[:,1])), axes_labels=["Xprime", "Vprime"], plotjoined=True)

#4.1.7
a=10
t=srange (0, 50, 0.1)
var ("X", "V")
sol= desolve_odeint((V, -X-(a*V^3-V)), ics= [5, 10], dvars=[X, V], times=t)
list_plot(list(zip(t, sol[:,0])), axes_labels=["time", "Vprime"]) + list_plot(list(zip(t, sol[:,1])), color="red")

#4.1.7
a=10
t=srange (0, 50, 0.1)
var ("X", "V")
sol= desolve_odeint((V, -X-(a*V^3-V)), ics= [5, 10], dvars=[X, V], times=t)
list_plot(list(zip(sol[:,0], sol[:,1])), axes_labels=["Xprime", "Vprime"], plotjoined=True)

#4.1.7
The parameter of a has an inverse relationship with the output of X' and V'. As you increase a, the subsequent time series and trajectories decrease. 

#4.2.1
@interact
def system (n=8):
    var("H", "P", "G")
    t=srange (0, 100, 0.1)
    k1=k2=k3=0.2
    Hprime= (1/(1+G^n))-(k1*H)
    Pprime= H-(k2*P)
    Gprime= P-(k3*G)
    sol= desolve_odeint([Hprime, Pprime, Gprime], ics=[1, 1, 1], dvars=[H, P, G], times=t)
    p= list_plot(list(zip(t, sol[:,0])), plotjoined=True, axes_labels= ["Time", "Hormone Secretion"]) + list_plot(list(zip(t, sol[:,1])), plotjoined=True, color="red") + list_plot(list(zip(t, sol[:,0])), plotjoined=True)
    show (p)

#4.2.2
var("H", "P", "G")
t=srange (0, 100, 0.1)
k1=k2=k3=0.2
n=8
Hprime= (1/(1+G^n))-(k1*H)
Pprime= H-(k2*P)
Gprime= P-(k3*G)
sol= desolve_odeint([Hprime, Pprime, Gprime], ics=[10, 20, 15], dvars= [H, P, G], times=t)
list_plot(sol)

#4.2.2
var("H", "P", "G")
t=srange (0, 100, 0.1)
k1=k2=k3=0.2
n=8
Hprime= (1/(1+G^n))-(k1*H)
Pprime= H-(k2*P)
Gprime= P-(k3*G)
sol= desolve_odeint([Hprime, Pprime, Gprime], ics=[15, 25, 20], dvars= [H, P, G], times=t)
list_plot(sol)

#4.2.2
var("H", "P", "G")
t=srange (0, 100, 0.1)
k1=k2=k3=0.2
n=8
Hprime= (1/(1+G^n))-(k1*H)
Pprime= H-(k2*P)
Gprime= P-(k3*G)
sol= desolve_odeint([Hprime, Pprime, Gprime], ics=[20, 30, 50], dvars= [H, P, G], times=t)
list_plot(sol)

#4.2.3
@interact
def system (n=(1, 10)):
    var("H", "G")
    t=srange (0, 100, 0.1)
    k1=k2=0.2
    Hprime= (1/(1+G^n))-(k1*n)
    Gprime= H- (k3*G)
    sol= desolve_odeint([Hprime, Gprime], ics=[1, 1], dvars= [H, G], times=t)
    p=list_plot(list(zip(t, sol[:, 0])), plotjoined=True, axes_labels=["Time", "Hormone Secretion"]) + list_plot(list(zip(t, sol[:, 1])), plotjoined=True, color= "red")
    show (p)