@sarra Your mutiple posts on the subject still do not answer my question about the requirement number four.

Let me put the question in a different way.

Suppose that someone tells you that the point Lt is at the location [3.5 4.2] relative to the orthogonal vectors A and B on the top-plate. Wouldn't you ask "But what are A and B"?

Whoever told you about the coordinates [3.5 4.2], should have also supplied you with the vectors A and B, otherwise that data is meaningless.

The first three specifications translate to the following Maple statements:

P := 13 * < 0, sin(10*Pi/180), cos(10*Pi/180) >;
N := < sin(18*Pi/180), 0, cos(18*Pi/180) >;
Lb := < 7, 5, 0 >;

The fourth specificaion is incomplete because it describes Lt relative to the vectors A and B, but A and B can be any pair of orthogonal vectors in the top plate.  Need to be more specific there.

@asceduardo The attached worksheet illustrates a very pedestrian approach to a multiple scale analysis of the differential equation that you have supplied.  It is incomplete since I have not applied the initial conditions, but you may complete it yourself if you feel sufficiently motivated.

multiscale.mw

Many thanks, Kitonum and vv, for your solutions.  I had expected pdsolve() to produce an answer in one step, but until that happens, the alternatives offered by you can be quite workable substitutes.

I have a few comments regarding your question.

1. The method of multiple scales would make sense when there is
a scale parameter, usually epsilon, in the equation or its
initial/boundary conditions.  Your equation x'' + x = 0 has
no epsilon in it, and you have not provided initial/boundary
conditions, therefore seeking a solution with epsilon is not
really meaningful.

2. In your proto-algorithm, you take
    t = T0 + e*T1 + e^2*T2
    x = e*X1(T1) + e^2*X2(T1) + e^3*X3(T1)
(here I am writing e for epsilon).

I don't think you mean that.  It's likely that you meant
    x = e*X1(T0,T1,T2) + e^2*X2(T0,T1,T2) + e^3*X3(T0,T1,T2)

3. It may be possible to do what you want through
PDEtools:-dchange(), but I haven't used that function and
I don't know how to help you use it.  If I were doing this
problem, I would do it in Maple exactly how I do it by
hand on paper, and that would not be very difficult.

In view of point #1 above, I cannot illustrate the steps,
but if you provide a more realistic problem, I may give
it a try.

To print a dot instead of an asterisk for multiplication:

@Carl Thank you very much for the detailed and clear explanation of `?[]`.

@Carl That's very nice.

Your use of the `?[]` operator is new to me.  I see that `?[]`(x,[a]) yields x[a].
How does that work and where is it documented?

@vv You are right; taking rij hat to be a unit vector in the rij direction, then both formuals are correct.

@Carl Actually the formula in the component form is correct.

It's the formula in the vector/bold form that needs to be fixed.  There, the r_ij in the denominator should be r_ij² to make the force proportional to inverse distance, as postulated.

If we let z be the expression of interest, then 1/expand(1/z) yields B immediately.  I don't know why expand(z) does not do the same thing.

@Axel Vogt The general solution in terms of hypergeom function is quite nice. Is it possible to simplify it when q is an integer?

Direct evaluations of the integral with arbitrary integers q (positive or negative) produce answers in terms of elementary functions.  I wonder if those may be obtained from the general formula.

@MapleUser2017 Perhaps someone else can jump in and show how to use the packages that you have indicated.

The first error message appears at the place where you do:

It would be a good idea to inspect what the expression for H_n(f) looks like before attempting to plot it.  When you do that, you will see that there are several undefined variables in it.  What is, for instance, N_p?

@acer Thanks for this very nice solution.  I would have never thought of it.