First, I would almost never use the trapezoidal rule to compute numerical integrals of smooth integrands when you can evaluate the integrand wherever you want. It is exponentially faster to use something like Gaussian quadrature. Moreover, packages like QuadGK can handle vector-valued integrands directly (rather than integrating each row separately):
using QuadGK
A = [1 2; 3 4]
b = [1;2]
t_final = 10
quadgk(t -> exp(A*t)*b, 0, t_final)[1]
gives
2-element Vector{Float64}:
3.7347743661398145e22
8.164742103848492e22
which is accurate to at least 12 digits.
However, in this particular case, the integral is exactly A^{-1} (e^{At} - I) b:
julia> using LinearAlgebra
julia> A \ ((exp(A*t_final) - I) * b)
2-element Vector{Float64}:
3.7347743661399164e22
8.164742103848739e22
(Of course, you can speed this up further as noted above by @Oscar_Smith by using a package to compute e^{At} b directly rather than computing e^{At} first.)