@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.

Hi Preben, that's rather curious.  I would have expected pdsolve to examine the PDE itself but as you have noted it appears that it merely looks at the boundary conditions!

Oscar, I still don't know what you mean by f.  Earlier you said it was the radial velocity, but now it seems that it's the pressure, or perhaps something else.

@adel-00  A phase portrait consists of a collection of orbits.  A system of differential equations has infinitely many orbits.  You don't want to sketch all of them, do you?

You need to specify a finite subset of the orbits that you wish to sketch.  You specify each orbit by an initial condition.  No one gives you the initial conditions—you are responsible for picking them.  The choice of initial conditions is very much problem-dependent.  You use your insight/understanding/taste to pick a good set.

Here is an example to get you started:
 

restart;
with(DEtools):
de1 := diff(x(t),t) = y(t);
de2 := diff(y(t),t) =  -x(t) - y(t);
de3 := diff(z(t),t) = -2*z(t);
ic := seq([x(0)= 1, y(0)=h, z(0)=1], h=-1..1, 0.25),
      seq([x(0)=-1, y(0)=h, z(0)=1], h=-1..1, 0.25);
DEplot3d([de1, de2, de3], [x(t), y(t), z(t)], t=0..10, [ic],
        linecolor=blue, thickness=1);

@oscmh3  I am trying to make sense of the equations.  Here are things that puzzle me.

1. You refer to f as "radial velocity".  To me that means velocity in the direction of a ray from the origin.  Is that what you mean?  In a tornado the velocity would be perpendicular to that direction.
2. If f is indeed radial velocity, then the Dirac(x-1) in the second initial condition says that the initial radial velocity at x=1 is very large, which would imply that particles located there will overtake their neighboring particles.  That's rather odd—in continuum mechanics one generally assumes that a material cannot penetrate itself.
3. In any case, the equation of motion derived from Newton's law would be expressed in terms of acceleration, which would be the first derivative of f with respect to y.  I don't understand why your equation involves the second derivative of f with respect to y.  The latter is the derivative of the acceleration!
4. You say that you want to model an explosion on a tornado's axis.  I am unable to connect that with the presence of both Dirac(0.1 - x) and Dirac(x - 1) in your equations.  What do the 0.1 and 1 signify?

@oscmh3 The eq2 in your original post looks suspicious.  For the second term in it you have:

I suspect it is supposed to be:
You may want to check that.