You might get more interest in this question on the Maxima mailing list. See: http://maxima.sourceforge.net/maximalist.html
For the first version of Lutus example 2, I get:
y(t) = m*sin(%omega*t)/(%omega^2*c^2*r^2+1)
-%omega*c*m*r*cos(%omega*t)/(%omega^2*c^2*r^2+1)
+(y(0)*%omega^2*c^3*r^3+%omega*c^2*m*r^2+y(0)*c*r)*%e^-(t/(c*r))
/(c*r*(%omega^2*c^2*r^2+1))$
which is the same as the expected result, if y(0) = 0. However, I don't see where that is assumed.
After atvalue(y(t),t=0,init_val)
, I get the same result as Lutus, namely:
y(t) = m*sin(%omega*t)/(%omega^2*c^2*r^2+1)
-%omega*c*m*r*cos(%omega*t)/(%omega^2*c^2*r^2+1)$
I am working with Maxima 5.31.1, built with Clisp, on Linux.