#
# Solve a spring system with 4 springs and 3 masses.
#
#
#                           Maple Midterm Project #2
#                  Harmonic Oscillators (Spring Systems)
#                                   by Mike Pilant
--------------------------------------------------------------------------------
# 
# Restart Maple...
# 
> restart;\

--------------------------------------------------------------------------------
# 
# Set parameters for four spring, three mass system:
# 
> N :=4;

                                     N := 4
# 
# Set spring stiffnesses here:
# 
> K[1] := 1.0:\
K[2] := 1.0:\
K[3] := 1.0:\
K[4] := 1.0:\

# 
# Set masses here:
# 
> M[1] := 1.0:\
M[2] := 2.0:\
M[3] := 3.0:
# 
# Set Forcing terms:
# 
> F[1] := t -> 0.0:\
F[2] := t -> 0.0:\
F[3] := t -> 0.0:\

# 
# Set Boundary conditions (corresponding to fixed walls):
# 
> x[0] := t -> 0.0:\
x[N] := t-> 0.0:\

--------------------------------------------------------------------------------
# 
# Set Initial conditions:
# 
> inits := x[1](0)=1.0, D(x[1])(0) = 0.0, \
             x[2](0) = 0.0, D(x[2])(0) = 0.0,\
             x[3](0) = 0.0, D(x[3])(0) = 0.0 :\

--------------------------------------------------------------------------------
# 
# Define governing equations (second order system):
# 
> inits := {}: eqs := {}: fcns := {}:\
for j from 1 to N-1 do\
  eq[j] := M[j]*diff( x[j](t), t$2 ) = K[j+1]*(x[j+1](t)-x[j](t)) - K[j]*(x[j](t)-x[j-1](t)) + F[j](t);\
  inits := inits union {x[j](0) = 0.0, D(x[j])(0) = 0.0};\
  eqs := eqs union {eq[j]};\
  fcns := fcns union {x[j](t)};\
od:\

# Now apply non-zero initial condition for one spring...
# 
> inits := inits union {x[1](0) = 1.0}:\

--------------------------------------------------------------------------------
# 
# Solve the system,
# 
> sol := dsolve( eqs union inits, fcns, method=laplace):\

--------------------------------------------------------------------------------
# 
# Define solution functions:
# 
> fns := {}:\
for k from 1 to N-1 do\
  z[k] := unapply( rhs(sol[k]), t);\
  fns := fns union {z[k]};\
od:\

--------------------------------------------------------------------------------
# 
# Plot solutions:
# 
> plot(fns);\