restart:
with(plots):
de := diff(y(t),t,t)+2*b*diff(y(t),t)+omega^2*y(t)=sin(alpha*t);
Undamped Case
b := 0;
de;
yt := A*sin(alpha*t)+B*cos(alpha*t);
subs(y(t)=yt, de);
expand(%);
collect(%, [cos(alpha*t), sin(alpha*t)]);
solve({-B*alpha^2+omega^2*B=0, -A*alpha^2+omega^2*A=1}, {A,B});
subs(%, yt);
de2 := subs(omega=1, alpha=Pi/5, de);
dsolve({de2, y(0)=-1, D(y)(0)=0}, y(t));
yy := unapply(rhs(%), t);
plot(yy(t), t=0..20*Pi, numpoints=200);
de3 := subs(alpha=1, omega=1, de);
yt := t*(A*cos(t)+B*sin(t));
subs(y(t)=yt, de3);
expand(%);
yy := subs(A=-1/2, B=0, yt);
dsolve({de3, y(0)=-1, D(y)(0)=2}, y(t));
rs := unapply(rhs(%), t);
plot(rs(t), t=0..20*Pi);
Damped Case
b := 'b';
de;
yt := A*sin(alpha*t)+B*cos(alpha*t);
subs(y(t)=yt, de);
expand(%);
collect(%, [sin(alpha*t), cos(alpha*t)]);
eqs := {-A*alpha^2-2*b*B*alpha+omega^2*A=1, -B*alpha^2+2*b*A*alpha+omega^2*B=0};
solve(eqs, {A,B});
amp := subs(%, sqrt(A^2+B^2));
simplify(%);
subs(omega=1, %);
amp := %;
plot3d(amp, b=0.1..2, alpha = 0..2,axes=boxed);
de4 := subs(b=0.2, omega=1, de);
dsolve({de4, y(0)=0, D(y)(0)=0}, y(t));
sol := rhs(%);
animate(plot, [sol, t=0..10*Pi], alpha = 0..2, frames=100);