human
I am working on the file "DifferentialEquations.ipynb".
The first task is:
Observations

We need to solve a system of equations on each time step.
The cost to solve the linear system is $O(n^3)$ for a system of $n$ equations.
The Jacobian matrix is $n\times n$, so requires $n^2$ storage when everything else requires $O(n)$.
I used finite differencing to compute the Jacobian, but many problems have efficient ways to compute an exact Jacobian (it just takes some programming).
Backward Euler is stable (as predicted by theory), but not very accurate.
What happens when you make the time step smaller or larger?

Midpoint Method
What if instead of evaluating the function at the end of the time step, we evaluated in the middle of the time step using the average of the endpoint values.  After all, something similar improved accuracy for our numerical integration...
$$ \tilde y(h) = y(0) + h f\left(\frac h 2, \frac{\tilde y(h) + y(0)}{2} \right) $$
Can you write Python code for it?
 Give me a Python code to compute the Jacobian, but many problems have efficient ways to compute an exact Jacobian (it just takes some programming)?
assistant
#Here is Python code to compute the Jacobian, but many problems have efficient ways to compute an exact Jacobian (it just takes some programming).

def ode_midpoint(f, y0, tfinal=1, h=0.1):
    y = y0.copy()
    t = 0
    thist = [t]
    yhist = [y0]
    while t < tfinal:
        h = min(h, tfinal - t)
        # Solve x = y + h f(x)
        def residual(x):
            return x - h * f(t+h/2, (x + y)/2) - y
        y = solve_newtonfd(residual, y)
        t += h
        thist.append(t)
        yhist.append(y.copy())
    return numpy.array(thist), numpy.array(yhist)

thist, yhist = ode_midpoint(test.f, y0, h=.2, tfinal=50)
pyplot.figure()
pyplot.plot(thist, yhist, '*')
pyplot.plot(thist, test.y(thist, y0))
pyplot.title('Midpoint')

y0 = numpy.array([.2])
pyplot.figure()
for test in tests:
    thist, yhist = ode_midpoint(test.f, y0, h=1, tfinal=20)
    pyplot.plot(thist, yhist, '*', label=repr(test)+' Midpoint')
    pyplot.plot(thist, test.y(thist, y0), label=repr(test)+' exact')
pyplot.legend(loc='upper right')
pyplot.title('Cosine')