Maple code for a rectangular pulse

m: mass, r: damping constant, k: spring constant,
a: constant in the impulse function ,
eps: width of a rectangular pulse , n: number of frames

> impulse_func:=proc(m,r,k,a,eps,n)
local h,dlign1,dlign2,sol1,sol2,plt1,plt2,pltd,plt,txt,txtd,kloss,e;
with(plots):
alias(u=Heaviside):
h:=eps->(u(t-a)-u(t-a-eps))/eps;
dlign1:=eps->m*diff(y(t),t,t)+r*diff(y(t),t)+k*y(t)=5*h(eps);
dlign2:=m*diff(y(t),t,t)+r*diff(y(t),t)+k*y(t)=5*Dirac(t-a);
print(dlign1(eps), epsilon=eps);
print(dlign2);
sol1:=eps->simplify(dsolve({dlign1(eps),y(0)=0,D(y)(0)=0},y(t)));
sol2:=dsolve({dlign2,y(0)=0,D(y)(0)=0},y(t));
e:= textplot([6,1.3,convert([101],bytes)], font=[SYMBOL,12]):
txt:=eps->textplot([6.9,1.3,cat(` = `,convert(evalf(eps,2),string))]);
txtd:=textplot([4,1.4,`Response to Dirac's delta function`],align={RIGHT,ABOVE});
pltd:=eps->plot(h(eps),t=0..5*Pi,color=green,thickness=2):
plt1:=eps->plot(rhs(sol1(eps)),t=0..5*Pi,color=red,thickness=2):
plt2:=plot(rhs(sol2),t=0..5*Pi,color=blue,thickness=2):
plt:=eps->display(pltd(eps),txt(eps),e,txtd,plt1(eps),plt2);
display(seq(plt(eps-(eps-0.1)/n*i),i=0..n),insequence=true,args[7..nargs]);
end: