import numpy as npfrom scipy.integrate import odeintimport matplotlib.pyplot as plt# function that returns dz/dtdef model(z,t,u):    x = z[0]    y = z[1]    dxdt = (-x + u)/2.0    dydt = (-y + x)/5.0    dzdt = [dxdt,dydt]    return dzdt# initial conditionz0 = [0,0]# number of time pointsn = 401# time pointst = np.linspace(0,40,n)# step inputu = np.zeros(n)# change to 2.0 at time = 5.0u[51:] = 2.0# store solutionx = np.empty_like(t)y = np.empty_like(t)# record initial conditionsx[0] = z0[0]y[0] = z0[1]# solve ODEfor i in range(1,n):    # span for next time step    tspan = [t[i-1],t[i]]    # solve for next step    z = odeint(model,z0,tspan,args=(u[i],))    # store solution for plotting    x[i] = z[1][0]    y[i] = z[1][1]    # next initial condition    z0 = z[1]# plot resultsplt.plot(t,u,'g:',label='u(t)')plt.plot(t,x,'b-',label='x(t)')plt.plot(t,y,'r--',label='y(t)')plt.ylabel('values')plt.xlabel('time')plt.legend(loc='best')plt.show()