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.

Assuming that column separators consist of exactly five spaces as stipulated—the OP's example violates this—then we may use StringSplit(), as in:

restart;
A:=Array():
k:=1:
fid:=fopen("/tmp/data", READ):
line := readline(fid):
while line <>0 do
          L := StringTools:-StringSplit(line, "     ");
          A(k):= [ L[1], parse(L[2]), parse(L[3]) ];
          line:=readline(fid);
          k:=k+1;
od:
fclose(fid);
M:= Matrix(k-1,3, (i,j) -> A[i][j] );

which yields:

I haven't examined the details of your calculations, but I am alarmed by any numeric calculation that involves numbers such as 10^27 and 10^(-11).

In solving differential equations numerically, we have to be acutely aware of discretization and floating point truncation errors that are naturally associated with such procedures.  It's very easy to obtain nonsensical results if you ignore these.

How do we avoid such pitfalls?  Rule number 1 is that you should use appropriate units.  The meter is a horribly bad unit for expressing distances between planets.  In modeling the Solar System, a good unit of length is the Astronomical Unit (AU) which is the average distance between the Earth and the Sun.  One AU is approximately 1.5x10^11 meters.  When you express the differential equation of motion the planets in terms of AU units, the distances will be smallish numbers which will be treated much more kindly under discretization and floating point truncations of the numerical solver.

While you are at it, you should also choose appropriate units for mass and time.  A kilogram is a good unit if you want to weigh a bag of potatoes.  It's a totally wrong unit if you want to express the mass of a planet.  For the time unit choose the length of a typical orbital period which would be of the order of a year, rather than a second.

Having selected the units of length, mass, and time, convert the gravitational constant G into those units.  If all is well, G will be a moderate-sized number rather than something like 6x10^(-11).

You should worry about Maple only after you are done with these preliminary preparations.

You don't need Maple for such a trivial problem.  Just solve the equation by hand and see what the solution looks like.

@acer After you have calculated A and B, you may obtain the plot immediately through:

plot([solve(ee,T), x, x=A..B]);

No need to call plottools:-reflect.

@Ronan ...let's see how his blog develops.

@Vaporwave thank you very much for your help.  The page you have found is exactly what I was looking for.  I wish MaplePrime's web tools worked more intelligently.

@Ronan It was a false alarm.  What I had was correct.  I have posted a new version which exapnds some of the comments imbedded in the worksheet.

Look up the help page for DEplot in the DEtools package.  It contains several examples.