Software 6.5.1 (Laplace transforms in SymPy). Transforms of standard functions, a shifted step, and the solution of Example 6.4.3 by partial fractions and inversion:

>>> import sympy as sp
>>> t, s = sp.symbols('t s', positive=True)
>>> a, b = sp.symbols('a b', positive=True)
>>> sp.laplace_transform(t**3, t, s, noconds=True)
6/s**4
>>> sp.laplace_transform(sp.exp(-a*t) * sp.cos(b*t), t, s, noconds=True)
(a + s)/(b**2 + (a + s)**2)
>>> sp.laplace_transform(sp.Heaviside(t - 2), t, s, noconds=True)
exp(-2*s)/s
>>> Y = 1 / (s * (s**2 + 3*s + 2))         # y'' + 3y' + 2y = 1 from rest
>>> sp.apart(Y, s)
1/(2*(s + 2)) - 1/(s + 1) + 1/(2*s)
>>> sp.inverse_laplace_transform(Y, s, t)
1/2 - exp(-t) + exp(-2*t)/2

Reading the output. The transforms agree with Theorem 6.4.2: \(3!/s^4\), the shifted cosine \((s + a)/((s + a)^2 + b^2)\), and the delayed step \(e^{-2s}/s\). apart gives the partial fractions of Example 6.4.3, and the inverse transform returns \(\frac 12 - e^{-t} + \frac 12e^{-2t}\). Declaring t and s positive avoids the conditions on \(\Real s\) that SymPy would otherwise carry.

MATLAB/Octave: syms t s; laplace(t^3) and ilaplace(1/(s*(s^2 + 3*s + 2))).

Software 6.5.2 (Characteristic roots and undetermined coefficients). A resonant case of Theorem 6.2.1, and an equation with a complex pair of characteristic roots:

>>> y = sp.Function('y'); u = sp.symbols('u', real=True)
>>> sp.dsolve(y(u).diff(u, 2) - 2*y(u).diff(u) + y(u) - sp.exp(u), y(u))
Eq(y(u), (C1 + u*(C2 + u/2))*exp(u))
>>> sp.dsolve(y(u).diff(u, 3) - y(u), y(u))
Eq(y(u), C3*exp(u) + (C1*sin(sqrt(3)*u/2) + C2*cos(sqrt(3)*u/2))*exp(-u/2))

Reading the output. For \(y'' - 2y' + y = e^u\), the forcing exponent \(1\) is a double root of \(p(\lambda ) = (\lambda - 1)^2\), so \(m = 2\) and the particular solution is \(\frac 12u^2e^u\), the term \(u^2/2\) inside the bracket. For \(y''' = y\), the roots of \(\lambda ^3 = 1\) are \(1\) and \(-\frac 12 \pm \frac {\sqrt 3}{2}i\), giving \(e^u\) and \(e^{-u/2}\) times \(\cos \) and \(\sin \) of \(\frac {\sqrt 3}{2}u\), as Theorem 6.1.2 predicts.

MATLAB/Octave: roots([1 0 0 -1]) for the characteristic roots.

Software 6.5.3 (The convolution formula). The oscillator of Figure 9 forced by a unit pulse on \([1, 3]\), solved by the convolution integral of Theorem 6.3.1 and by direct numerical integration:

>>> import numpy as np
>>> from scipy.integrate import solve_ivp, quad
>>> w0, z = 2.0, 0.1                       # y'' + 2 z w0 y' + w0^2 y = g(t)
>>> wd = w0 * np.sqrt(1 - z**2)
>>> G = lambda t: np.exp(-z * w0 * t) * np.sin(wd * t) / wd * (t >= 0)
>>> g = lambda t: 1.0 * (1 <= t <= 3)     # a rectangular pulse
>>> conv = lambda t: quad(lambda s: G(t - s) * g(s), 1, min(t, 3))[0]
>>> rhs = lambda t, x: [x[1], g(t) - 2*z*w0*x[1] - w0**2*x[0]]
>>> num = solve_ivp(rhs, (0, 8), [0, 0], max_step=0.01, rtol=1e-10,
...                 atol=1e-12, dense_output=True)
>>> for T in (2.0, 3.0, 5.0, 8.0):
...     print(T, round(conv(T), 5), round(num.sol(T)[0], 5))
...
2.0 0.31452 0.31452
3.0 0.37458 0.37458
5.0 -0.12393 -0.12393
8.0 -0.10297 -0.10297

Reading the output. The two methods agree to five decimals at every time. The convolution needs only the impulse response and a quadrature, and it gives the response to any forcing once \(G\) is known; this is how linear systems are characterised experimentally, by measuring their impulse or step response. Note the response after the pulse ends at \(t = 3\): the oscillator keeps ringing, with amplitude about \(0.12\) at \(t = 5\), decaying like \(e^{-0.2t}\).

MATLAB/Octave: impulse(tf(1, [1 0.4 4])) and lsim with the Control System Toolbox.