My PhD was on the topic of delay differential equations (DDEs) and, at the time, I wrote lots of utilities to do things like calculate the stability of linearised equations. All that code has been lost to the mists of time but I still have call for it every now and again (e.g., for real-time dynamic substructuring where I need to model the delays in the control loop).

Here is some code in Julia that calculates the stability of a DDE in the form

$$ x’(t) = Ax(t) + Bx(t-\tau) $$

where $A$ and $B$ can be (square) matrices, i.e., systems of DDEs. If you’ve got a nonlinear model, you just need to linearise around an equilibrium first.

The algorithm implemented is that of Pseudospectral Differencing Methods for Characteristic Roots of Delay Differential Equations; D. Breda, S. Maset, and R. Vermiglio, SIAM Journal on Scientific Computing 2005 27:2, 482-495. It has the property that it approximates the characteristic roots of the DDE closest to the origin best.

using LinearAlgebra

"""
    cheb_mesh(T::Type, N::Integer)
    cheb_mesh(N::Integer)

Return the Chebyshev nodes of degree `N` on the interval [0, 1].
"""
cheb_mesh(T::Type, N::Integer) = (-cos.((T(π) * (0:N)) ./ N) .+ 1) ./ 2
cheb_mesh(N::Integer) = cheb_mesh(Float64, N)

"""
    cheb(T::Type, N::Integer)
    cheb(N::Integer)

Return a tuple `(x, D)` where `x` contains the Chebyshev nodes of degree `N` on the interval
[0, 1] using the numerical type `T` and `D` is the Chebyshev differentiation matrix on the
same interval.

## Example - differentiate a sinusoid

    julia> (x, D) = cheb(10) ;

    julia> D * sinpi.(x) ≈ π * cospi.(x)
    true
"""
function cheb(T::Type, N::Integer)
    if N == 0
        return (ones(T, 1), zeros(T, 1, 1))
    else
        x = cheb_mesh(T, N)
        c = [2; ones(T, N - 1, 1); 2] .* (-1) .^ (0:N)
        dx = x .- x'
        D = (c * (1 ./ c)') ./ (dx + I)
        D[diagind(D)] .-= vec(sum(D; dims=2))
        return (x, D)
    end
end
cheb(N::Integer) = cheb(Float64, N)

"""
    dde_stability(A, B, τ; N=10)

Calculate (a finite number of) the characteristic roots of the linear delay differential
equation

    x'(t) = Ax(t) + Bx(t-τ)

where `A` and `B` can be (square) matrices. `N` controls the order of the Chebyshev
polynomial approximation to the infinitesimal generator; higher values of `N` result in a
more accurate approximation. A property of the algorithm used is that characteristic roots
closest to the origin are approximated best.

## Reference

Pseudospectral Differencing Methods for Characteristic Roots of Delay Differential
Equations; D. Breda, S. Maset, and R. Vermiglio, SIAM Journal on Scientific Computing 2005
27:2, 482-495
"""
function dde_stability(A, B, τ; N=10)
    @assert size(A) == size(B)
    @assert size(A, 1) == size(A, 2)
    T = eltype(A) == Int ? Float64 : eltype(A)
    n = size(A, 1)
    (x, D) = cheb(T, N)
    M = kron(-D, Matrix(I / convert(T, τ), n, n))
    M[1:n, 1:n] .= A
    M[1:n, (n + 1):(end - n)] .= 0
    M[1:n, (end - n + 1):end] .= B
    return eigen(M)
end

# Tests
# 1. x'(t) = 2x(t) - exp(1)*x(t - 1) has a double root at λ = 1
# 2. x'(t) = -x(t - τ) goes unstable at τ=π/2
# 3. x'(t) = α*A*x(t) + (1-α)*A*x(t-τ)
#    A = [
#          50  284   41   23   50   32
#        -280  -46  -19  -37  -10  -28
#          35   -1   26  143   35   17
#          5   -31 -139  -22    5  -13
#          20  -16   11   -7 -115  137
#         -10  -46  -19  -37 -145 -163
#    ]
#    with α=0.9 is stable for 3π/(270*(2*0.9 - 1)) ≤ τ ≤ 4π/270
#    that is 0.04363323129985824 ≤ τ ≤ 0.046542113386515455