The finite element method described above is most successful in treating boundary value problems.  When addressing initial value problems, while the finite element method could be applied (and sometimes is), other techniques are frequently preferable.  FlexPDE uses a variable-order implicit backward difference method (BDM) as introduced by C.W. Gear.  In most cases, second order gives the best tradeoff between stability, smoothness and speed, and this is the default configuration for FlexPDE.  This method fits a quadratic in time to each nodal value, using two known values and one future (unknown) value.  It then solves the coupled equations for the array of nodal values at the new time.  By looking backward one additional step, it is possible to infer the size of the cubic term in a four-point expansion of the time behavior of each nodal value.  If these cubic contributions are large, the timestep is reduced, and if extreme, the current step repeated.