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

It's likely that the fourth boundary condition is misstated because it's redundant—it's implied by the thrid condition.  Check!

Additionally, you need to specify the domain over which the PDE is to be solved.

Furthermore, why do you expect this PDE to have a solution at all?  What is the motivation behind it?

@bronze45629 The subjet in this correspondence has become a moving target, making the communication difficult.  So let me give the complete calculations here to avoid any confusion:

restart;
with(DEtools):
satODE1 :=
    diff(vx(t),t) = -G*M*x(t)/(x(t)^2+y(t)^2)^(3/2),
    diff(vy(t),t) = -G*M*y(t)/(x(t)^2+y(t)^2)^(3/2),
    diff(x(t),t) = vx(t),
    diff(y(t),t) = vy(t);
ics := [x(0)=1, y(0)=0, vx(0)=0, vy(0)=1];
subsys := subs({G=1,M=1},[satODE1]);
DEplot(subsys, [x(t),y(t),vx(t),vy(t)], t=-2..2,
    [ics], scene=[x(t),y(t)], scaling=constrained);

See if that helps.

As to your question regarding brackets, I don't want to start writing a mini-lesson here on Maple's syntax.  Just look over the code I have provided here and modify it a little at a time as needed.

@bronze45629 Suppose you want to plot not one, but several orbits at once.  In that case you would specify;

ics := [ x(0)=a1, y(0)=b1, vx(0)=c1, vy(0)=d1],
          [ x(0)=a2, y(0)=b2, vx(0)=c2, vy(0)=d2],
          [etc];

You see that the square brackets group the initial conditions of indivitual oribits.

When you call DEplot, you would do:

DEplot(...., [ics], ....);

The brackets here group the multiple initial conditions into a single list.

In your case you are plotting only a single orbit, therefore you would do;

ics := [ x(0)=a1, y(0)=b1, vx(0)=c1, vy(0)=d1];
DEplot(...., [ics], ....);

Do you see that?

@bronze45629 In your original post you had:

ics:=[x(0)=1, y(0)=0,vx(0)=0,vy(0)=1];

In your second post you have removed the brackets.  Restore them.

Perhaps you mean

   plot([seq([u[i], p[i]], i = 1 .. N)]);

So far as I know, Maple is not equipped to handle elliptic equations.  But it may be able to handle time-dependent ones.  Try the time-dependent case to see what happens.

@Christopher2222 In case of irregular spacing, tomleslie's version will still work.