import numpy as np
import matplotlib.pyplot as plt
def f(i,t,y) :
	k=np.zeros([2,1])
	k[0] = y[0]*(2.0-0.01*y[1])
	k[1] = y[1]*(0.01*y[0]-1)
	return k[i]

def F(yold,x,t,h) :
	z=np.zeros([2,1])
	z[0] = x[0] - yold[0] - h/2*f(0,t-h/2,x[:])
	z[1] = x[1] - yold[1] - h/2*f(1,t-h/2,x[:])
	return z[:]

def J(t,h,y) :
	z=np.zeros([2,2])
	z[0,0] = 1 - h/2*(2-0.01*y[1])
	z[0,1] = h/2*0.01*y[0]
	z[1,0] = -h/2*0.01*y[1]
	z[1,1] = 1 - h/2*(0.01*y[0]-1)
	return z[:,:]

def newton(t,yold,h,d) :
	w=np.zeros([2,1])
	yeuler=np.zeros([2,1])
	w[:] = yold[:]
	for k in range(3) :		
		Jinv = np.zeros([2,2])
		Jinv = np.linalg.inv(J(t,h,w[:]))
		yeuler = w - np.dot(Jinv,F(yold[:],w[:],t,h))
		w[:] = yeuler[:] 
	return yeuler[:]


def myeuler(a,b,y0,M,d) :
	h=float((b-a))/M
	ynew=np.zeros([2,1])
	yexp=np.zeros([2,1])
	yold=np.zeros([2,1])
	yoldexp=np.zeros([2,1])
	data1 = np.zeros([2,M+1])
	data1[:,0] = y0
	data2 = np.zeros([2,M+1])
	data2[:,0] = y0
	time = np.zeros([M+1])
	time[0] = a
	t=a
	yold[:,0] = y0
	yoldexp[:,0] = y0
	for n in range(M) :
		yeuler = newton(t+h,yold[:],h,d)
		for i in range(d) :
			ynew[i] =2*yeuler[i] - yold[i]
			yexp[i] = yoldexp[i] + h*f(i,t,yoldexp[:])
			data1[i,n+1] = ynew[i]
			data2[i,n+1] = yexp[i]
		t = t + h
		time[n+1] = t
		yold[:] = ynew[:]
		yoldexp[:] = yexp[:]	
	return data1[:,:], data2[:,:], time[:]
a=0.0
b=20.0
t0=a
d=2
y0=np.array([300, 150])
M=1000
(data1, data2, time) = myeuler(a,b,y0,M,d)
fig = plt.figure()
plt.plot(time[:],data2[0,:],'r-',label='y1-exp euler')
plt.plot(time[:],data1[0,:],'b--',label='y1-middle')
plt.title('Explicit Euler vs Middle Point Method')
plt.legend()
plt.show()

fig = plt.figure()
plt.plot(time[:],data2[1,:],'r-',label='y2-exp euler')
plt.plot(time[:],data1[1,:],'b--',label='y2-middle')
plt.title('Explicit Euler vs Middle Point Method')
plt.legend()
plt.show()

fig = plt.figure()
plt.plot(time[:],data1[0,:],'r--',label='prey')
plt.plot(time[:],data1[1,:],'b--',label='predator')
plt.legend()
plt.show()

fig = plt.figure()
plt.plot(data2[0,:],data2[1,:],'r--',label='predator population (Exp Euler)')
plt.plot(data1[0,:],data1[1,:],'b--',label='predator population (Midpoint Rule)')
plt.legend()
plt.show()