--Copyright The Numerical Algorithms Group Limited 1994.
-- Input file to solve linear second order differential equations
-- with driving force A*cos(wt).

-- The equation is solved at y(0)=y'(0)=0 for 
--    (i)   an overdamped, forced motion example  c**2-4*k*m > 0
--    (ii)  critically damped                     c**2-4*k*m = 0
--    (iii) underdamped                           c**2-4*k*m < 0

-- The resulting equations are then plotted.

)clear all
y := operator y
deq := m*D(y x, x, 2) + c*D(y x, x) + k*(y x) = A * cos (w * x)

solve(deq, y, x=0, [0,0])                     -- takes a few minutes
ex:=%
eval(%, [c=6,k=5,m=1,A=6*sqrt(5),w=sqrt(5)])  -- c**2-4*k*m = 16
draw(%,x=0..15,title=="Overdamping")

eval(ex, [k=5,m=1,A=6*sqrt(5),w=sqrt(5)])     -- c**2-4*k*m = 0
limit(%,c=sqrt(20),"right")                   -- otherwise division by zero
draw(%,x=0..15,title=="Critically Damped")

eval(ex, [c=2,k=5,m=1,A=6*sqrt(5),w=sqrt(5)]) -- c**2-4*k*m = -16
trigs %
rule1 := rule sin(-x) == - sin(x)
rule2 := rule cos(-x) ==   cos(x)
rule2 rule1 %% (-3)
draw(%,x=0..15,title=="Underdamped")
