Symbolic Math Toolbox 

Consider an ideal massspringdamper system with mass m (in kg), spring constant k (in N/m) and damping coefficient R (in Ns/m).
This system can be described with the following formula: . We use the ode function to define the equation of motion and the initial conditions.
f1 := ode({m*x''(t) + R*x'(t) + k*x(t), x(0) = 0 , x'(0) = 1}, x(t))
Assuming m = 10, R = 1 and k = 10, the equation can be rewritten as. We substitute these values in our differential equation expression and solve the equation.
f2 := subs(f1, m = 10 , R = 1 , k = 10)
f3 := solve(f2);
As we can see above, the solution has an exponential decay term and a sinusoidal term. This implies that the system will show an underdamped response to a perturbation, as can be seen in the Figure 1 below.
p := plot::Point2d([t, f3[1]], t=0..50):
c := plot::Curve2d([t,f3[1]], t = 0..a, a= 0..50):
plot(p,c);
Note that Figure 1 was customized by modifying line properties and adding a title within the notebook interface’s object browser. These properties could have also been specified programmatically in the plot commands above.
Figure 2 shows an animation of the oscillating mass. This animation was created in a separate notebook and copied into this notebook. The ability to copy graphics and animations between notebooks is useful when users wish to document supporting analysis without having the code displayed in the notebook.
The code used to create this animation is available here.
We will now analyze the frequency response of our massspringdamper system. We begin by deriving the Laplace transform of our massspringdamper system, and then generate Bode plots (magnitude vs. frequency and phase vs. frequency). We perform these tasks within a custom function called FreqResp, which outputs the bode plots for a given ODE input. The Bode plot indicates that the cutoff frequency of our massspringdamper system is approximately 1 rad/sec.
fn := 10*x''(t) + x'(t) + 10*x(t)  h(t):
FreqResp(fn)