I need to plot the function f(k) which is defined as a square root of the time integral of the square of the Euclidean norm of the system state vector x(t) during the first 10 seconds after the start of movement. Input is matrix A, vector b and inital state x0 for a system of differential equations dx/dt = Ax + bu. u is called control signal and is of the form u = kx. Dimension of vector x is 2. If system is unstable f(k) should return -1.

Here is what I’ve tried (without plotting, just wanted to get correct function f(k)):

import numpy as np
from scipy.integrate import solve_ivp
from scipy.integrate import quad
from numpy.linalg import norm


def f(A, b, k, x0):
    M = A + b * k
    eigenvalues = np.linalg.eig(M)[0]
    if eigenvalues[0] > 0 and eigenvalues[1] > 0:
        return -1
    else:
        def odefun(t, x):
            x = x.reshape([2, 1])     
            dx = M @ x      
            return dx.reshape(-1)

        # solve system
        x0 = x0.reshape(-1)
        sol = solve_ivp(odefun, [0, 10], x0)['y']
        print(sol.shape) # (2, 41)
        return quad(norm(sol), 0, 10)


def contra(A, b, x0):
    vals = []
    for k in range(10):
        vals.append(f(A, b, k, x0))
    return vals

A = np.array(([1, 2], [3, 4]))
b = np.array([5, 6])
x0 = np.array([7, 8])
print(contra(A, b, x0))

But I’m getting this error where I call function quad:

ValueError: invalid callable given

Could somebody please give a hint of how to write this correctly? Thanks in advance.