%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

# s = state
# t = time
def derivs(s,t):
m = s[0]
c = s[1]
dmdt = +0.1*m - 0.15*c*m
dcdt = -.1*c + 0.15*c*m
return np.array([dmdt,dcdt])
def EulerStep(s, t, derivs, dt):
"""
EulerStep is a generic implimentation of the Euler Algorithm.
Notice that it knows *nothing* about the state, nor the physics
of the problem (which is handled by the 'derivs' function
that is passed in by the programmer).
"""
return s + derivs(s,t)*dt
def HeunStep(s, t, derivs, dt):
"""
HeunStep is a generic implimentation of the Heun Method.
Notice that it knows *nothing* about the state, nor the physics
of the problem (which is handled by the 'derivs' function
that is passed in by the programmer).
"""
# print(s)
# print(t)
f1=derivs(s,t)
f2=derivs(s+f1*dt,t+dt)
return s + 0.5*(f1+f2)*dt
s = np.array([10,10]) ## mice, cats
dt = 0.01 # pick a reasonable time step
t = 0.0 # start at t=0.0
tf = 10.0 # stop in 500 sec
c = 0
m = 0
tlist = [t]
clist = [c]
mlist = [m]
s = HeunStep(s,t,derivs,dt)
t = t + dt
while t <= tf:
s = HeunStep(s, t, derivs, dt)
t += dt
tlist.append(t)
mlist.append(s[0])
clist.append(s[1])

first_half_tlist = tlist[0:(len(tlist)//2)] # first half of tlist
second_half_tlist = tlist[(len(tlist)//2):-1] # second half of tlist
c_heun_first = np.array(clist[0:(len(tlist)//2)])
c_heun_second = np.array(clist[(len(tlist)//2):-1])
m_heun_first = np.array(mlist[0:(len(tlist)//2)])
m_heun_second = np.array(mlist[(len(tlist)//2):-1])
plt.plot(first_half_tlist,c_heun_first,'r-',label="Cats")
plt.plot(first_half_tlist,m_heun_first,'b-',label="Mice") # first half
plt.legend()
plt.title("Cats and mice relationship system")
plt.xlabel("time (s)")
plt.ylabel("Number of cats and mice")
plt.show()

plt.plot(second_half_tlist,c_heun_second,'r-',label="Cats")
plt.plot(second_half_tlist,m_heun_second,'b-',label="Mice") # first half
plt.legend()
plt.title("Cats and mice relationship system")
plt.xlabel("time (s)")
plt.ylabel("Number of cats and mice")
plt.show()