I think that by increasing the time scales to three or higher instead of two, your results will be a little better, But it will be more complex solutions. Also 'ode45' method is not suitable for solving eigenvalue problems and only solving at the exact eigenvalue can be trusted.

Thank you for your answer.

I agree with you that increasing the time scales may result in higher accuracy in the time domain response.

The ode45 in the mentioned paper however, was not used for solving the eigenvalue problem. It was used to obtain the nonlinear frequency response from the solution of the time-domain problem. 

There are several error sources in the use of the Method of Multiple Time Scales (MMTS).

1. The expansion, like MOST expansion methods does not provide a convergent series. Usually, the series is an asymptotic series.  it does not converge.  Rather, it provides a NUMERICAL approximation to the solution.  That approximation is limited in validity to some time (in most cases, only to times of order 1/epsilon).  Hence, if one insists on using the results of MMTS for times that are too long, the results are bound to deviate from the exact solution.

2.  Updating the fundamental frequency.  If this is done through some order in epsilon (epsilon^(n+1)) then an error will accumulate of order epsilon^n beyond times of order 1/epsilon.

3. Truncating the expansion in some order.  The more orders one adds, the better is the accuracy (still, limited in time!).

4. If one goes beyond first order, one MUST check for the compatibility of the solvability conditions in different orders.  For example, if one goes as far as order epsilon^2 in the expansion, one must check for the compatibility of the first- and second-order solvability conditions.  If one does not fulfill the constraint that this consistency check imposes on the solution, the results may be nonsensical.

This requires that, in the differential equation in a given order, one must include the general solution of the homogeneous equation.  Often, the consistency between the solvability condition in that order and the solvability condition in the next order imposes a constraint on the general solution of the homogeneous equation.

Now to your paper.  I notice that you expand s(t) through second order, while you expand all other entities only through first order.  Also, you use only T0 and T1.  This is inconsistent.  Everything has to be expanded to the same order: q(t), r(t) and s(t), and, concurrently, you must use T0, T1, T2 and T3.  Then, going as far  as O(epsilon^2) in the solution, bd as far as O(epsilo^3) in the fundamental frequency, you MUST include in q1, r1 and a1 free terms that represent the general solution of the homogeneous equation, and find what the requirement for consistency between the first- and second-order solvability conditions imposes on these free terms.

Thank you Professor Zarmi for your comprehensive answer.

The reason for having s(t) in one order higher than the other two variables is to have it present in the the same equation of suppression for the system. Since s(t) it is in one order lower, if it the same order of epsilon is considered, it will not show up in the equations. 

I assume moving to a higher order of epsilon and including all terms will help a more precise response as you mentioned.

Dear Ehsan,

Having gone through your paper again, I insist on my previous statement.  For example, in your Eq. (17), you couple s(t) and q(t).  Hence BOTH have to be expanded through the same order in epsilon for consistency.  If you want to work through O(epsilon^2), you must expand both s(t) and q(t) through the same order.  If you want the O(epsilon^2) approximation to be valid through times of O(1/epsilon), you must expand everything through O(epsilon^3), so as to find the updating of the fundamental frequency through O(epsilon^3) due to its T3 dependence.

At the end, you will have the solution expanded through O(epsilon^2), and the fundamental frequency - through O(epsilon^3).

This is the general "rule of thumb" for using MMTS.  You expand everything through O(epsilon^(n+1)), so as to find the dependence of the fundamental frequency on T0, T1, T2, ....., T(n+1).  Then you sum the result for the solution through O(epsilon^n) (one order LESS than in the fundamental frequency, which is to be found through O(epsilon^(n+1)).  This yields an approximation with an error of O(epsilon^(n+1)) for times of O(1/epsilon).

Yair

Dear Prof. Zarmi

Thank you for your complete and comprehensive answer.

I greatly appreciate it. 

The best way to view the method of multiple scales in oscillation problems is to look at it from the view point of a time varying coordinate transform.  That way one can avoid the physical indeterminacy of the multiple times and multiple amplitudes.  See my online book at http://www.maths.adelaide.edu.au/anthony.roberts/modelling.php : click on the 3rd topic  " The normal form of oscillations give their amplitude and frequency".

There are many advantages:  by rephrasing as a coordinate transform you can view the coordinate transform as exact and the simple dynamics in the new variables as an approximation that you can tightly control.  Or you can view the coordinate transform as approximately linking the simple system and the original.  Or you can view the coordinate transform together with the simple system as being exactly for some system close to the original. 

Furthermore, this coordinate transform view empowers you to break free from the tyranny of having to scale all variables and parameters to a single epsilon.  If you wish you can approximate systems with multiple small parameters with the smallness being independent of each other.  Or you can link their smallness with one epsilon.  The choice is yours.

BTW for those interested in the method of multiple scales in space-time problems instead of just time, I have recently proposed analysis to rigorously bound the error in the modelling of the evolution of slow space structures.  See "Macroscale, slowly varying, models emerge from the microscale dynamics in long thin domains" at http://arxiv.org/abs/1310.1541

Ehsan, if you'd like you can send me in gmail details of the problem (if this one is not an essentially nonlinear, or non-smooth, from which let me be free, - beginning from a generating solution), it's better, some working version (type of Maple-protocol or similar - much more fine than any text) for understanding. I'll try your eqs explicitly.