Hi Anurag, Jason
I think using ODE module isn't compulsory in for finding solution of the equation x' = Ax + Bu , I've tried it by using sympy.solvers.ode.systems.linodesolve but it was giving weird result like this , I don't understand what is Dummy terms here.

I've another idea to find the solution of this differential equation by using Inverse Laplace Transform and by using it the process becomes more simpler. This approach will be more simple and I think inverse_laplace can be used in many places inside of control.
The final solution in this case for a Forced response with input will be
x(t) = inv_laplace([SI - A]**-1)X(0) + inv_laplace([SI-A]**-1 * B * U(S))
X(0) is the initial response it can be also denoted as X(t_0) here.
Here we just need to find the inverse of SI - A once and put it where ever needed.