restart;
with(plots); with(plottools);
omega := 1;
alpha := 1.5;
b := 0.1;
y0 := 0;
yp0 := 0.1;
de := diff(y(t),t,t)+2*b*diff(y(t),t)+omega^2*y(t)=sin(alpha*t);
ss := dsolve({de, y(0)=y0, D(y)(0)=yp0}, y(t));
yy := unapply(rhs(ss), t);
tmax := 8*Pi;
yscale := 2;
numsteps := 50;
dt := tmax/numsteps;
for i from 0 to numsteps do tt:= i*dt; p[i] := display([line([-1,-yscale*y0],[-1,yscale*y0], color=green), point([-1, yy(tt)], symbol=circle, color=red), plot(yy(t),t=-0.01..tt,color=red), line([-1, yy(tt)],[tt, yy(tt)], color=blue)]); od:
display([seq(p[i], i=0..numsteps)], view=[-1..tmax, -yscale*abs(y0)..yscale*abs(y0)],insequence=true);
de2 := diff(y(t),t,t)+omega^2*y(t)=0;
ss2 := dsolve({de2, y(0)=y0, D(y)(0)=yp0}, y(t));
yy2 := unapply(rhs(ss2),t);
display([plot(yy(t), t=0..6*Pi, color=red), plot(yy2(t), t=0..6*Pi, color=green)]);