NDSolve - backward difference method

There is always a mess in various names of various numerical methods. Mathematica's NDSolve is a highly sophisticated black box. It is quite difficult to uncover the method which is actually used during computation.
However, here I have example solutions of one equation using default and "Backward Difference Formula" of order 1 (i.e. Implicit Euler) Method. As we will see later, the "default" method in this case is "Backward Difference" of higher order.

eqn=-u0 \[Omega]^2 Sin[t \[Omega]] + (g x[t])/r + 2 b x'[t] + x''[t];
vals = {b -> 0.001, r -> 0.42578125, g -> 9.81, u0 -> 0.01}; 

ns1 = NDSolve[{0 == eqn /. vals /. \[Omega] -> 3.14, 
    x[0] == 0.001, x'[0] == 0.001}, x, {t, 0, 200}][[1]];

ns2 = NDSolve[{0 == eqn /. vals /. \[Omega] -> 3.14, 
    x[0] == 0.001, x'[0] == 0.001}, x, {t, 0, 200}, 
   Method -> {"BDF", "MaxDifferenceOrder" -> 1}, MaxSteps -> 10^6, 
   PrecisionGoal -> 5, AccuracyGoal -> Automatic][[1]]

Plot[{x[t] /. ns1, x[t] /. ns2}, {t, 100, 120}, AspectRatio -> 1/5, 
  ImageSize -> 800, PlotStyle -> {Green,Red}]

Here is the result: green is the default solution, red is the 1st order BDF.

I have a similar plot, where the blue and red curves come from numerical solutions using gsl routines implementing so called Gear's stepping method: gsl_odeiv_step_gear1 and gsl_odeiv_step_gear2 (it seems, that the current version of gsl has abandoned these routines, it uses gsl_odeiv2_step_msbdf instead, see docs). Blue is the second order solution (together with points of accepted solution), red is the fisrt order solution (there was too many small steps to show them). The green line of the "default" Mathematica solution line hides behind the blue "gear2" curve.

StepMonitor and InterpolationFunction

One could assume that InterpolatingFunction resulting from a numerical algorithm like NDSolve would have the same node distribution as the algorithm.
And it seems to be true.

{s, steps} = 
  Reap[NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, 
    y, {x, 0, .05}, StepMonitor :> Sow[{x, y[x]}][[1]]]];

Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
coords = First[InterpolatingFunctionCoordinates[y /. s[[1]]]]
steps[[1, All, 1]]

The both coords and steps give the same values:

{5.13196*10^-6, 0.0000102639, 0.00339028, 0.00677029, 0.0101503, 
0.0151503, 0.0201503, 0.0251503, 0.0301503, 0.0351503, 0.0401503, 
0.0450752, 0.05}

However, the Mesh option of Plot shows another set of points (green):

stpl = ListPlot[steps, 
   PlotStyle -> Directive[PointSize[Large], Red]];
Show[
 Plot[Evaluate[y[x] /. s], {x, 0, .016}, Mesh -> All, 
  MeshStyle -> Directive[PointSize[Medium], Green]], stpl]

And it seems that the Mesh directive has no connection with the structure of InterpolatingFunction. Let's make larger Plot and cut of a part which corresponds to the previous figure:

stpl = ListPlot[steps, 
   PlotStyle -> Directive[PointSize[Large], Red]];
Show[
 Plot[Evaluate[y[x] /. s], {x, 0, .05}, Mesh -> All, 
  MeshStyle -> Directive[PointSize[Large], Green]], stpl, 
 PlotRange -> {{0, 0.016}, {1, 1.0085}}]

Plot axes labels, legends, fonts etc

Insert text into plot using Epilog command - here y comes to the vertical axis, \omega to horizontal axis and some legend to the right upper corner of the figure:

Epilog -> {
  Text["y", Scaled[{0.05, .95}], Background -> White,
                                 BaseStyle -> {FontSize -> 11, Bold}],
  Text["\[Omega]", Scaled[{0.95, .1}], 
                                 Background -> White, 
                                 BaseStyle -> {FontSize -> 11, Bold}],
  Text["b=" <> ToString[b] <> ",u=" <> ToString[u], 
                        Scaled[.98 {1, 1}], {1, 1}, Background -> White]
 }

Alternatively, use the Text's align feature, but this does not work well at the y axis:

Epilog -> {
  Text[" \[Xi] ", Scaled[{0.025, 1}], {-1, 1.2}, 
   Background -> White],
  Text[" \[Omega] ", Scaled[{1, 0}], 1.2 {1, -1}, Background -> White,
    BaseStyle -> {FontSize -> 14}],
  Text["\[Beta]=" <> 
    ToString[b], Scaled[.98 {1, 1}], {1, 1}, Background -> White]
  }

Setting up the default style for texts in Plots:

bst = {FontSize -> 12, 
            FontFamily -> "Helvetica", 
            FontWeight -> Bold, 
            FontSlant -> Italic};
SetOptions[Plot,BaseStyle->bst];
SetOptions[ListPlot,BaseStyle->bst];
SetOptions[Text, BaseStyle -> bst];

Align axes in several plots: It is necessary to manually allocate space for axes labels via , ImagePadding -> {{left, right}, {bottom, top}}
 

{
GraphicsColumn[{Plot[Sin[x], {x, 0, 2 \[Pi]}], 
                Plot[10^-3 Sin[x], {x, 0, 2 \[Pi]}]}],
ip = {{40, 0}, {5, 5}};
GraphicsColumn[{Plot[Sin[x], {x, 0, 2 \[Pi]}, ImagePadding -> ip], 
                Plot[10^-3 Sin[x], {x, 0, 2 \[Pi]}, ImagePadding -> ip]}]
}