Software 1.7.1 (Symbolic solution with SymPy). dsolve recognises the standard types; checkodesol substitutes the answer back into the equation. For the exact equation of Example 1.4.4 the potential is built by the two integrations in the proof of Theorem 1.4.2:

>>> import sympy as sp
>>> t = sp.symbols('t'); y = sp.Function('y')
>>> lin = sp.dsolve(y(t).diff(t) + 2*y(t) - sp.exp(-t), y(t)); lin
Eq(y(t), (C1*exp(-t) + 1)*exp(-t))
>>> sp.checkodesol(y(t).diff(t) + 2*y(t) - sp.exp(-t), lin)
(True, 0)
>>> ber = sp.dsolve(y(t).diff(t) - y(t) + y(t)**3, y(t)); ber[1]   # y > 0
Eq(y(t), sqrt(-exp(2*t)/(C1 - exp(2*t))))
>>> x, v = sp.symbols('x v')               # exact: M dx + N dv = 0
>>> M, N = 2*x*v, x**2 + 3*v**2
>>> sp.diff(M, v) == sp.diff(N, x)          # the exactness test
True
>>> F = sp.integrate(M, x)                  # M in x, then the v part
>>> F = F + sp.integrate(sp.simplify(N - sp.diff(F, v)), v); F
v**3 + v*x**2

Reading the output. The linear equation \(y' + 2y = e^{-t}\) has general solution \(y = e^{-t} + C_1e^{-2t}\), printed in the factored form \((C_1e^{-t} + 1)e^{-t}\), and checkodesol confirms it with residual \(0\). For the Bernoulli equation \(y' = y - y^3\), the positive branch is \(y = \bigl (e^{2t}/(e^{2t} - C_1)\bigr )^{1/2}\), which the substitution \(w = y^{-2}\) of Theorem 1.5.1 produces by hand; every positive solution tends to \(1\). The potential of the exact equation is \(x^2y + y^3\), as in Example 1.4.4 (with \(v\) for \(y\)). Asking dsolve for the exact equation itself returns the three roots of the cubic \(y^3 + x^2y - C = 0\) in Cardano’s form, a reminder that the implicit form \(F(x, y) = C\) is often the better answer.

MATLAB/Octave: syms y(t); dsolve(diff(y) + 2*y == exp(-t)) with the Symbolic Math Toolbox (or the symbolic package in Octave).

Software 1.7.2 (Numerical solution with solve_ivp). The logistic equation with \(r = 1\), \(K = 100\), \(y_0 = 10\), solved numerically and compared with the exact solution of Example 1.3.3:

>>> import numpy as np
>>> from scipy.integrate import solve_ivp
>>> r, K, y0 = 1.0, 100.0, 10.0
>>> exact = lambda t: K / (1 + (K / y0 - 1) * np.exp(-r * t))
>>> sol = solve_ivp(lambda t, y: r * y * (1 - y / K), (0, 10), [y0],
...                 rtol=1e-8, atol=1e-10, dense_output=True)
>>> ts = np.array([0, 2, 4, 6, 10])
>>> np.round(np.c_[ts, sol.sol(ts)[0], exact(ts)], 4)
array([[ 0.    , 10.    , 10.    ],
       [ 2.    , 45.0853, 45.0853],
       [ 4.    , 85.8486, 85.8486],
       [ 6.    , 97.8178, 97.8178],
       [10.    , 99.9592, 99.9592]])
>>> sol.nfev, np.max(np.abs(sol.sol(ts)[0] - exact(ts)))
(314, 8.665153856668439e-07)

Reading the output. The numerical and exact values agree to four decimals: the largest difference at the five times is \(8.7 \times 10^{-7}\), within the requested tolerances. The solver, an adaptive Runge–Kutta method of order \(5(4)\) that is described in the last section, used \(314\) evaluations of \(f\). The population passes \(K/2 = 50\) just after \(t = 2\), at \(t^* = \log 9 = 2.20\), and reaches \(99.96\) by \(t = 10\). dense_output=True returns an interpolant, sol.sol, so the solution can be evaluated at any time, not just at the solver’s own steps.

MATLAB/Octave: opts = odeset(’RelTol’, 1e-8); then [t, y] = ode45(@(t, y) y.*(1 - y/100), [0 10], 10, opts).

Software 1.7.3 (Equilibria and the harvesting bifurcation). The equilibria of \(y' = y(1 - y) - h\) and the signs of \(f'\) at them, for three harvest rates:

>>> u, h = sp.symbols('u h', real=True)
>>> f = u * (1 - u) - h
>>> eq = sp.solve(f, u); eq
[1/2 - sqrt(1 - 4*h)/2, sqrt(1 - 4*h)/2 + 1/2]
>>> for hv in [0.16, 0.25, 0.3]:
...     roots = [complex(e.subs(h, hv)) for e in eq]
...     slopes = [complex(sp.diff(f, u).subs({u: rt, h: hv})) for rt in roots]
...     real = [abs(rt.imag) < 1e-12 for rt in roots]
...     print(hv, [round(rt.real, 3) for rt, ok in zip(roots, real) if ok],
...           [round(s.real, 3) for s, ok in zip(slopes, real) if ok])
...
0.16 [0.2, 0.8] [0.6, -0.6]
0.25 [0.5, 0.5] [0.0, 0.0]
0.3 [] []

Reading the output. For \(h = 0.16\) the equilibria are \(0.2\) and \(0.8\), with \(f' = 0.6\) (unstable) and \(-0.6\) (stable), as in Figure 1. At \(h = 0.25\) they coincide at \(0.5\), where \(f' = 0\) and the linear test of Theorem 1.6.1(iv) fails. For \(h = 0.3\) the roots are complex and there are no equilibria at all: the population declines to zero from any starting value.

MATLAB/Octave: roots([-1 1 -0.16]) for the equilibria at \(h = 0.16\).