def euler_step(f,t0,t1,y):

"""

Compute next value for Euler's method ODE solver

Args:

f: right-hand side function that gives rate of change of y

y: float value of dependent variable

t0: float initial value of independent variable

t1: float final value of independent variable

Returns:

1d array of dependentant variables updated after the timestep

"""

h = t1 - t0

return y + h*f(t0, y)

def rk2_step(f,t0,t1,y):

"""

compute next value for Euler's method ODE solver

Args:

f: right-hand side function that gives rate of change of y

y: float value of dependent variable

t0: float initial value of independent variable

t1: float final value of independent variable

Returns:

1d array of dependentant variables updated after the timestep

"""

h = t1 - t0

#compute euler estimate for half step

y1 = y + 0.5*h*f(t0, y)

t1 = t0 + 0.5*h

#compute midstep rate estimate and take full step using midstep rate estimate

return y + h*f(t1, y1)

def rk4_step(f,t0,t1,y):

"""

compute next value for 4th-order Runge-Kutta method ODE solver

Args:

f: right-hand side function that gives rate of change of y

y: float value of dependent variable

t0: float initial value of independent variable

t1: float final value of independent variable

Returns:

1d array of dependentant variables updated after the timestep

"""

pass

def rk_solve(f,t,y0,order):

"""

Runge-Kutta solver for systems of 1st order ODEs

Args:

f: right-hand side function that gives rate of change of y

y0: numpy array of initial float values of dependent variable

t: numpy array of float values of independent variable

order: int order of RK solver with allowed values [1,2,4]

Returns:

2D numpy array of float values of dependent variable at each time step

"""

if order == 1:

method = euler_step

elif order == 2:

method = rk2_step

elif order == 4:

method = rk4_step

else:

raise ValueError("rk_solve `order` must be 1, 2, or 4")

y = [y0]

for i in range(t.size - 1):

y.append(method(f, t[i], t[i + 1], y[i]))

return np.array(y)

def compute_solution(f, tmin, tmax, y0, steps, order=2):

'''

Set up time intervals and initial conditions and compute numerical solution

Args:

f: function describing rates of change

eps: float parameter in rate function

Returns:

t: 1D numpy float array of independent variable values

y: 2D numpy float array of dependent variable values

'''

t = np.linspace(tmin, tmax, steps + 1)

y = rk_solve(f,t,y0,order)

return t, y

def simulate_RKF45(f, tmin, tmax, y0, steps):

'''

Function to set up proper arguments to perform simulation using ivp_solve.

See scipy documentation for details.

Note that ivp_solve may store solution data transposed from what we coded up previously.

Return:

t: 1D numpy array of independent variable values

y: 2D numpy array of dependent variable values

'''

pass

def writen_response_p3():

#3a Based on the plot of your numerical solution, describe the long-term steady-state behavior of the system.

#3b What happens when you try this simulation with 500 steps with a Van der Pol mu of 0.1? What is the steady-state behavior with 5000 steps?

#3c Compare and contrast the results of your RK4 and numpy's RK45 results. How many steps were used to get that agreement?

return """ANSWER: """