解答例

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

def func(x,a,b,c,d):#関数funcが入力xに対してその微分を返す。ただし入力x、出力dxは要素数2のベクトルとしている
    dx=np.zeros(2)
    dx[0]=a*x[0]-b*x[0]*x[1]
    dx[1]=c*x[0]*x[1]-d*x[1]
    return dx

a=1.0           
b=1             
c=1
d=1
dt=0.001           #時間刻み
t_max=30         #tの最大値

times=np.arange(0,t_max,dt)      #timesは出力される時刻を要素としたベクトル
total_steps=len(times)                 #total_stepsは総ステップ数

x=np.zeros((total_steps,2))           #（total_steps行、2列）の出力用の2次元行列であるxを用意
x[0,0]=2                                     #xの初期値を最初の要素に代入。x[0,0]とx[0,1]の最初の要素は時刻、次の要素は変数のインデックス
x[0,1]=1

for i in np.arange(total_steps-1):                 #オイラー法のメインループ
    x[i+1,:]=x[i,:]+func(x[i,:],a,b,c,d)*dt         #オイラー法の計算。ここでx[i,:]は、i番目の時刻の状態xを表す要素数2のベクトルになる

plt.plot(times,x[:,0],label="x_0") 
plt.plot(times,x[:,1],label="x_1")
plt.legend()
plt.xlabel("time")        
plt.ylabel("x")     
plt.show()