function sys!(dx,x,p,t)
mul!(dx, A, x)
mul!(dx, B, u, 1, 1)
nothing
end
this is the non-allocating version (depending a bit on where u is coming from).
If u is constant or piecewise constant, you can solve this differential equation by computing a single matrix exponential of exp([A B; 0 0]) and iterate the corresponding system
\begin{align}
x(k+1) = A_d x(k) + B_d u(k)
\end{align}