Refactor solve_ivp example in numerical integration module for clarity and completeness

Author: Daniel Precioso

modules/numerical-integration/index.qmd

@@ -28,52 +28,38 @@ Where:
 
 # SciPy's `solve_ivp`
 
-The `scipy.integrate.solve_ivp` function is the standard tool for solving ODEs in Python:
+The `scipy.integrate.solve_ivp` function is the standard tool for solving ODEs in Python. Try the following code:
 
 ```python
+import numpy as np
 from scipy.integrate import solve_ivp
 
-def dydt(t, y):
-    # Define your ODE system
-    return # derivative
-    
-sol = solve_ivp(dydt, t_span=(0, 10), y0=[initial_condition])
-```
-
-## Key Parameters
-
-- `fun`: The function defining the ODE system
-- `t_span`: Tuple `(t_start, t_end)` for integration interval
-- `y0`: Initial conditions (array)
-- `method`: Integration method (default: 'RK45')
-- `t_eval`: Specific time points to return solution
+def exponential_decay(t, y):
+     return -0.5 * y
 
-# Example: Exponential Growth
+t_span = [0, 10]
+y0 = [2, 4, 8]
+sol = solve_ivp(
+    fun=exponential_decay,
+    t_span=t_span, 
+    y0=y0)
 
-$$\frac{dy}{dt} = ky, \quad y(0) = y_0$$
+print(sol.t)
 
-```python
-import numpy as np
-import matplotlib.pyplot as plt
-from scipy.integrate import solve_ivp
-
-def exponential_growth(t, y, k):
-    return k * y
+print(sol.y)
+```
 
-k = 0.5
-y0 = [1.0]
-t_span = (0, 10)
-t_eval = np.linspace(0, 10, 100)
+**How does `solve_ivp` work?** Let's understand its parameters (copied from the [documentation](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html)):
 
-sol = solve_ivp(exponential_growth, t_span, y0,
-                args=(k,), t_eval=t_eval)
+- `fun`: Right-hand side of the system: the time derivative of the state $y$ at time $t$. The calling signature is `fun(t, y)`, where `t` is a scalar and `y` is an ndarray with `len(y) = len(y0)`. Additional arguments need to be passed if `args` is used (see documentation of `args` argument). `fun` must return an array of the same shape as `y`.
+   In our example, `exponential_decay` defines the ODE $\frac{dy}{dt} = -0.5y$, which models exponential decay $y(t) = y_0 e^{-0.5t}$.
+- `t_span`: Interval of integration $(t_0, t_f)$. The solver starts with $t=t_0$ and integrates until it reaches $t=t_f$. Both $t_0$ and $t_f$ must be floats or values interpretable by the float conversion function.
+- `y0`: Initial state. For problems in the complex domain, pass `y0` with a complex data type (even if the initial value is purely real).
 
-plt.plot(sol.t, sol.y[0])
-plt.xlabel('Time')
-plt.ylabel('y(t)')
-plt.title('Exponential Growth')
-plt.show()
-```
+You can also specify additional parameters:
+- `method`: Integration method to use. Common choices include `'RK45'` (default), `'RK23'`, `'DOP853'`, `'Radau'`, `'BDF'`, and `'LSODA'`.
+- `t_eval`: Times at which to store the computed solution, must be sorted and lie within `t_span`. If None (default), use points selected by the solver.
+- `args`: Additional arguments to pass to the user-defined functions. If, for example, `fun` has the signature `fun(t, y, a, b, c)`, then `args=(a, b, c)`.
 
 # Systems of ODEs