Alternative to @achresz's answer : change tools. Running :
X0, Y0, t, a0, b0, A0, a1, b1, A1, = SR.var('X0, Y0, t, a0, b0, A0, a1, b1, A1')
X, Y = function("X, Y")
de1 = X(t).diff(t) == a0*X(t)+b0*Y(t)-A0
de2 = Y(t).diff(t) == a1*X(t)+b1*Y(t)-A1
import sympy
%time foo = sympy.dsolve(*map(sympy.sympify, ([de1, de2], [X(t), Y(t)])))
%time tfoo=[u.rewrite("sin") for u in foo]
%time sfoo=[u._sage_() for u in foo]
%time stfoo=[u._sage_() for u in tfoo]
gives :
CPU times: user 1.78 s, sys: 8.09 ms, total: 1.79 s
Wall time: 1.79 s
CPU times: user 24.9 ms, sys: 56 µs, total: 24.9 ms
Wall time: 24.9 ms
CPU times: user 10 ms, sys: 0 ns, total: 10 ms
Wall time: 10 ms
CPU times: user 7.47 ms, sys: 3.89 ms, total: 11.4 ms
Wall time: 11.4 ms
The expressions are somewhat bulky, and neither Sympy nor Sage can lighten them significatively.
In this case, "the competition" (the gratis-but-not-free Wolfram Engine) does better. Running :
M1=mathematica("Sys={X'[t]==a0*X[t]+b0*Y[t]-A0, Y'[t]==a0*X[t]+b1*Y[t]-A1}")
M2=mathematica("TSol=Timing[FullSimplify[DSolve[Sys, {X[t], Y[t]}, t]]]")
gives
CPU times: user 215 µs, sys: 13 µs, total: 228 µs
Wall time: 12.2 s
(half a magnitude slower), but :
sage: M2[2]
{{X[t] -> (A1*b0 - A0*b1)/(a0*b0 - a0*b1) + E^(((a0 + b1)*t)/2)*
(C[1]*Cosh[(Sqrt[a0^2 + 4*a0*b0 - 2*a0*b1 + b1^2]*t)/2] +
(((a0 - b1)*C[1] + 2*b0*C[2])*
Sinh[(Sqrt[a0^2 + 4*a0*b0 - 2*a0*b1 + b1^2]*t)/2])/
Sqrt[a0^2 + 4*a0*b0 - 2*a0*b1 + b1^2]),
Y[t] -> (A0 - A1)/(b0 - b1) + E^(((a0 + b1)*t)/2)*
(C[2]*Cosh[(Sqrt[a0^2 + 4*a0*b0 - 2*a0*b1 + b1^2]*t)/2] +
((2*a0*C[1] - a0*C[2] + b1*C[2])*
Sinh[(Sqrt[a0^2 + 4*a0*b0 - 2*a0*b1 + b1^2]*t)/2])/
Sqrt[a0^2 + 4*a0*b0 - 2*a0*b1 + b1^2])}}
waaaay lighter...
Checking the identity of the solutions symbolically is problematic : (M2[2][1][1][2]-stfoo[0].rhs()._mathematica_()).FullSimplify()
doesn't return after several minutes.
Checking them numerically is (lazily) left as an exercise to the reader.
EDIT : Fixed a typo (mis-pasting) in Mathematica's code.
HTH,