Erratic outputs from a second order system when input is too small I'm modeling a simple Mass-Spring-Damper system to represent the torsional behavior of a micromirror. With references to some papers (this one mainly), I've constructed the model of a torsional mass spring damper system. 
Governing equation:
$$m\dfrac{d^2 \theta}{dt^2} + c \dfrac{d \theta}{dt} + k~\theta = T$$ 
Where T = Torque.
Calculated values (micromirror):
$m = 3.61*10^{-4} kg*um^2$ - -($3.61*10^{-16}) kg*m^2$
$c = 4.62*10^{-3} uN*um*s$ - - ($4.62*10^{-15}) N*m*s$
$k = 1.03 uN*um$ - - ($1.03*10^{-12}) uN*um$

The main issue I've been having is that below a certain input torque (~10^-6), my output seems to become unstable. I have shown examples of the good case (input is large enough), and the bad case (input is too small).
Input = ~10^-5 (in microns)

Input = ~10^-9 (in microns)

I'm stumped - I've been looking at this for a while without any solid grasp as to what is causing this to happen. I can increase the input by as much as I want and I will still see stable behavior, it is only at low input that this occurs. Any ideas what I could be doing wrong or neglecting?
 A: Your problem is poorly scaled: I don't know enough about ODE45 to try to debug what's going on under the hood, but any time you are passing in data to a black-box numerical method with quantities on the order of $10^{-21}$ (your torque of $10^{-9}$ micron microNewtons, in SI) you are asking for trouble.
You can try rescaling your problem so that your are working in gram-microns, and see if your issue goes away. Alternative, your black-box package is probably way overkill for such a simple ODE; here for instance is a snippet of code I wrote that solves it using Velocity Verlet:
#include <iostream>

void step(double &theta, double &dtheta, double dt, double m, double c, double k, double T)
{
    theta += dt*dtheta;
    dtheta += (T-c*dtheta-k*theta)*dt/m;
}

int main()
{
    // initial conditions
    double theta = 0;
    double dtheta = 0;

    // physical parameters
    double m = 3.61e-16;
    double c = 4.62e-15;
    double k = 1.03e-12;
    double T = 1e-21;

    // simulation parameters
    double dt = 0.001;
    double endt = 5;

    double curt = 0;
    for(int i=0; curt <= endt; i++)
    {
        curt += dt;
        double curT = 0;
        if(curt > 1.0)
            curT = T;
        step(theta, dtheta, dt, m, c, k, curT);
        std::cout << curt << " " << theta << std::endl;
    }
}

It works even without rescaling. The results I get for your torques are:

These are not identical to your results (right order of magnitude, at least) but I trust mine more ;) I leave it to you to double-check the constants, play with the timestep dt, etc.
