restart;with(plots):with(plottools):spring:=proc(x1,x2,y1,y2,n)
local kink,st,w;
st:=(y2-y1)/(n+1);
w:=(x2-x1)/2;
kink:=i->[w+x1,y1+i*st-st/2],[x2,y1+i*st-st/4],
[x1,y1+i*st+st/4],[w+x1,y1+i*st+st/2];
plot([[x1+w,y1],[x1+w,y1+st/2],seq(kink(i),i=1..n),
[x1+w,y2-st/2],[x1+w,y2]],
color=black,thickness=2);
end:motion:=proc(m,c,k,y0,v0,T,N)
local spring1,mass,pist,equil,cont,fluid,
mark,gplot,p2,frame,
y,dash,piston,springfun,mn,A,phi,
x1,x2,y1,y2,xscl,yscl,step;
step:=T/N;
if c^2>(k/m)
then y:=t->(y0*(c+sqrt(c^2-k/m))+v0)/2/sqrt(c^2-k/m)
*exp((-c+sqrt(c^2-k/m))*t)
+(y0*(-c+sqrt(c^2-k/m))-v0)/2/sqrt(c^2-k/m)
*exp((-c-sqrt(c^2-k/m))*t);
elif c^2=(k/m)
then y:=t->((v0+c*y0)*t+y0)*exp(-c*t);
else y:=t->exp(-c*t)*(y0*cos(sqrt(k/m-c^2)*t)
+(v0+c*y0)/sqrt(k/m-c^2)*sin(sqrt(k/m-c^2)*t));
fi;
mn:=evalf(minimize(y(t),[t=0..T]));
A:=evalf(maximize(y(t),[t=0..T])-mn);
mn:=min(mn,0);
xscl:=T/5;
yscl:=A/5;
spring1:=y1->spring(T+xscl/2,T+3*xscl/2,y1+yscl/2,A+yscl,5);
mass:=y1->rectangle([T+2*xscl/3,y1+yscl/2],[T+4*xscl/3,y1-yscl/2],color=red);
cont:=plot([[T+xscl/2,mn-2/3*yscl],[T+xscl/2,mn-A-7/3*yscl],
[T+3/2*xscl,mn-A-7/3*yscl],[T+3/2*xscl,mn-2/3*yscl]],
thickness=2,color=black);
fluid:=rectangle([T+xscl/2,mn-yscl],[T+3*xscl/2,mn-A-7/3*yscl],color=blue);
pist:=y1->plot([[T+xscl/2+.05*xscl,y1-A-3/2*yscl],
[T+3*xscl/2-.05*xscl,y1-A-3/2*yscl],
[T+xscl,y1-A-3/2*yscl],
[T+xscl,y1-yscl/2]],color=brown,thickness=3);
mark:=t->rectangle([t-.05*xscl,y(t)+.1*yscl],[t+.05*xscl,y(t)-.1*yscl],
color=red);
equil:=plot([[0,0],[T+2*xscl,0]],color=grey,thickness=1);
gplot:=plot(y(t),t=0..T,color=blue);
p2:=i->spring1(y(i*step));
if c>0
then frame:=i->display(gplot,equil,mark(i*step),mass(y(i*step)),
p2(i),fluid,cont,pist(y(i*step)),
view=[0..(T+1.75*xscl),
(mn-A-7/3*yscl-.05)..(A+yscl)],
axes=boxed);
else frame:=i->display(gplot,equil,mark(i*step),mass(y(i*step)),
p2(i),
view=[0..(T+1.75*xscl),-3/2*A..(A+yscl)],
axes=boxed);
fi;
display(seq(frame(j),j=0..N),insequence=true);
end:Here's the prototype for the motion procedure: motion(mass, damping, spring const, initial displacement, initial velocity, max. time, frames);mass=1, damping = .6, spring cosntant = 4motion(1,.6,4,2,0,8,30);mass = 1, damping = 2, spring constant = 4motion(1,2,4,2,0,5,30);mass = 1, damping = 4, spring constant = 4motion(1,4,4,2,0,5,30);No dampingmotion(1,0,4,2,0,8,30);