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.

@9009134 The attached worksheet shows how to solve Poisson's boundary value problem with finite differences.  It is not exactly the solution of your problem, but it is close enough so that you may use it as a template and modify it into yours.

I have commented some of the steps but I can't tell whether that would be comprehensible to you since I don't know about your mathematical background.  You may find details of finite differences in most numerical analysis textbooks, for instance, in Kincaid and Chaney.

Worksheet: FDM2.mw

Edit: Replaced the original FDM.mw by FDM2.mw.  The calculations are the same but I have added a few more comments.

@Carl Love and Kitonum, thank your very much for your answers to my question.  They are both very satisfactory.

@Kitonum That elimination method works but it's not ideal because requiring the twice differentiability of f1(t) is too much.

Those differential equations make sense even if f1(t) and f2(t) are discontinuous.  Solving the system through the Laplace transform requires only the integrability of f1 and r2.