Create a tangent approximation of a non-linear function f(x) about point a
    using a block lower triangular solver

    0 = f(x) = f(a) + J*x   # taylor series about a (if f(x) linear in x, then globally valid)
    J*x = -f(a)             # solve for x
    x = -J^{-1}f(a)         # but inverse is slow, so we use solve
    where J = df/dx

def tangent_approx(f: SYM, x: SYM, a: SYM = None, assert_linear: bool = False) -> Dict[str, SYM]:
    """
    Create a tangent approximation of a non-linear function f(x) about point a
    using a block lower triangular solver

    0 = f(x) = f(a) + J*x   # taylor series about a (if f(x) linear in x, then globally valid)
    J*x = -f(a)             # solve for x
    x = -J^{-1}f(a)         # but inverse is slow, so we use solve
    where J = df/dx
    """
    # find f(a)
    if a is None:
        a = ca.DM.zeros(x.numel(), 1)
    f_a = ca.substitute(f, x, a)  # f(a)
    J = ca.jacobian(f, x)
    if assert_linear and ca.depends_on(J, x):
        raise AssertionError('not linear')
    # solve is smart enough to to convert to blt if necessary
    return ca.solve(J, -f_a)