Here is a version of your method that uses solve instead of bfallroots.
Of course, this code will fail when solve fails to find the roots. When p has
symbolic roots (for example p = s^2 - (a+b)*s + a*b)), it might be
fun to detect this and try to return a conditional.
(%i11) L(p,f, s, t) := block([b,i : 0],
p : ratexpand(p),
b : solve(p,s),
b : xreduce('append, (map(lambda([s,k], makelist(s,i,1,k)), b, multiplicities))),
f : f / coeff(p, s, length(b)),
for bk in b do (
q : exp(rhs(bk) * t),
f : q * (integrate(f / q, t) + concat('c, i : i + 1))),
f)$
Try a double root:
(%i18) (s-5)^2 * (s-7)$
(%i19) ratexpand(%);
(%o19) s^3-17*s^2+95*s-175
Symbolic forcing term:
(%i20) L(%,g(t),s,t);
(%o20) %e^(5*t)*(integrate(integrate(%e^(2*t)*integrate(%e^(-7*t)*g(t),t),t),t)+(c1*%e^(2*t))/4+c2*t+c3)
And check:
(%i21) ratsimp(diff(%,t,3) - 17 *diff(%,t,2) + 95 * diff(%,t) - 175 * %);
(%o21) g(t)
--Barton
-----maxima-bounces at math.utexas.edu wrote: -----
>???(%i1)?L(p,f,m):=(fpprec:m,?b:bfallroots(p),?n:length(b),
>f:f/coeff(p,s,n),for?k?thru?n?do
>(q:exp(rat(rhs(b[k]))*t),f:trigrat(q*(integrate(f/q,t)+c[k]))),
>f:expand(bfloat(f)))$