You may want to have a look at DEplot3d.
The command you gave in your question should work if you replace DEplot with DEplot3d.
With the revised command you get a trajectory in phase space.

EOM is not an equation;  it is not of the form a=b;
If the intention was EOM = 0 then just change it to that.

You never defined ST, so how can you expect the select command to work?

Reading the Programming Guide Chaper 6 under Functional Operators: Mapping Notation we find:

Internally, a procedure created using operator notation is the same as any other procedure, except that it will have options operator, arrow.

To create a procedure where scoping rules are an issue it seems to me that you would have to enter it without using the arrow, i.e. as
p:=proc(...) option operator, arrow;  ... end proc;

Do you have an example in mind?

You should upload code which we can copy into Maple. Nobody is going to type this himself.

Taking a glance at the code there seems to be an implicit assumtion that w = (b-a)/h is an integer.
Since b is set at 10*n this means that (10*n-a)/h is supposed to be an integer.
If that is the case then changing the beginning to
for i from 1 to w do
....
would define x[0].
Try the following, where only the (implicitly defined) table x is defined.
The output is the sorted list of the table indices:
 

p:=proc(a,n,h) local b,w,x,i;
  b:=10*n;
  w:=(b-a)/h;
  x[w]:=10*n;
  for i from 1 to w do
    x[w-i]:=x[w-i+1]-h
  end do;
  sort([indices(x,nolist)]); #Getting the indices of the table x
end proc;

The exponential function exp, should be used like this exp(x), NOT like this e^(x).
I corrected that below.
You assign pi to 3.143. But you also use Pi. Below I have assumed that pi was intended as a (bad) approximation to Pi, so I have used Pi in all places.
So we have:
 

restart;
Digits:=50:
G := 6.6743*10^(-8);
a := 1.9501*10^24;
b := .3306;
c := 2.99792458*10^10;
d := 2.035;
#pi := 3.143; #Pi also appears below. I changed pi to Pi.
eps := 3.8220*10^35;
## g(r) and j(r) are not used in the following:
#g(r) = 1-s(r)/0.06123;
#j(r) = exp(-(1/2)*w(r))*(1-2*G*v(r)/(r*c^2))^(1/2); #Must use exp( ) not e^( )
sys1 := diff(v(r), r) = 4*Pi *r^2*eps/c^2, v(0)=0;
res1:=dsolve({sys1});
sys2:=diff(u(r), r) = -G*(eps+u(r))*(v(r)+4*Pi*r^3*u(r)/c^2)/(c^2*(r^2-2*G*r*v(r)/c^2)),u(0)=1.3668*10^34;
sys3:=diff(w(r), r) = 1.485232054*10^(-28)*(v(r)+4.450600224*10^(-21)*Pi *r^3*u(r))/(r^2-2*G*r*v(r)/c^2),w(0)=0;
## Now solve numerically
res2:=dsolve(eval({sys2,sys3},res1),numeric);
plots:-odeplot(res2,[[r,u(r)],[r,10^35*w(r)]],0..700000);
res2(688240);

The result of the last command is that w(r) = 0.0704472 ... and not as you "want"  w(688240)=-2.05684.
You cannot impose both conditions w(0)=0 and w(688240)=-2.05684.
 

You have a variety of choices. Here I give 3 methods:

restart; 
with(LinearAlgebra): 
A := Matrix(2,2,[-4,-2,3,1]); 
X:=<x1(t),x2(t)>;
ode:=diff~(X,t)=A.X; #The ode 
##Finding eigenvalues and corresponding eigenvectors (the columns in P): 
Lambda,P:=Eigenvectors(A); ## The general solution is
X = add(c[i]*exp(Lambda[i]*t)*Column(P,i),i=1..2); #FIRST
## That (in another form) could have been found by dsolve: 
dsolve(ode); 
X=subs(%,X); #SECOND
## You also have available MatrixExponential:
E:=MatrixExponential(A,t);
## The general solution in a third form: 
X = E.<c1,c2>; #THIRD

If you are not using a recent version of Maple, then dsolve won't take vector input.
In that case, replace ode by

ode1:=Equate(op(ode));
## and then as above:
dsolve(ode1);
X=subs(%,X);

Since it appears that you want to compare the numerical solution with the exact one you can try this:

restart;
Digits:=15;
ode := diff(w(x), x$6)+diff(w(x), x$4)+diff(w(x), x, x)-w(x) = -90;
bcs:=w(0) = 0, w(1) = 0, (D@@2)(w)(0) = 0, (D@@2)(w)(1) = 0, D(w)(0) = 0, D(w)(1) = 0; 
res:= dsolve({ode,bcs}, numeric);
plots:-odeplot(res,[x,w(x)]);
pol:=DEtools[de2diffop](lhs(ode),w(x),[R,x]); #The characteristic polynomial
rts:=solve(pol=0,R); #Roots not "nice"
evalf(rts);
## An impatient attempt at finding the exact solution:
sol:=timelimit(30,dsolve({ode,bcs})); #30 seconds expires
## I gave up on that also because the solution would take up a lot of space.
####
solG:=dsolve(ode): #The general solution
## Now imposing the boundary conditions:
eq1:=eval(rhs(solG),x=1):
eq2:=eval(diff(rhs(solG),x),x=1):
eq3:=eval(diff(rhs(solG),x,x),x=1):
eq4:=eval(rhs(solG),x=0):
eq5:=eval(diff(rhs(solG),x),x=0):
eq6:=eval(diff(rhs(solG),x,x),x=0):
## Use fsolve; it is fast:
S:=fsolve({eq1,eq2,eq3,eq4,eq5,eq6});
S:=simplify(fnormal(S,10));
sol:=simplify(evalc(evalf(eval(solG,S))));
plot(rhs(sol),x=0..1);
plots:-odeplot(res,[x,w(x)-rhs(sol)]); #Comparison

The reason is that there is no solution. You may also try some linear algebra:

M,b:=LinearAlgebra:-GenerateMatrix([A1, A2, A3, A4, A5, A6],[C1, C2, C3, C4, C5, C6]);
LinearAlgebra:-LinearSolve(M,b); #Inconsistent system, so no solution
## Details:
LinearAlgebra:-GaussianElimination(<M|b>);
### QUESTION:
What does the second part have to do with the first?
Consider
ode:=diff(w(x), x, x, x, x, x, x)+diff(w(x), x, x, x, x)+diff(w(x), x, x)+(1-2)*w(x) = -90;
dsolve(ode);
evalf(%);
## So why are you looking at
w(x)= C1*sinh(x)+C2*cosh(x)+C3*sin(x)+C4*cos(x)+C5*sin(x)+C6*cos(x)+90;
?
Notice that regardless of that question, the appearance of sin(x) and cos(x) twice each is strange!

Begin by deleting maxfun as the error message says.
Then you could try using approxsoln=[psi(u)=...] as an option.

My own experiments indicate that a solution looks like this:

Use the view option in Maple 2015 (which you are using, it is not needed in Maple 2016):

 plot3d(2*x/(x^2+y^2), x = -10 .. 10, y = -10 .. 10,view=-5..5);

For commenting out several lines use   (*  ....   *)  as in

(*

line 1
line 2
line 3
*)
 

You have a system of two odes that is third order in f and first order in theta.

Since 3+1 = 4, you need 4 boundary conditions.

In your other version, which I deleted, the problem is the same.

In the help page in Maple 2016, dsolve/details, you find under Optional Arguments:

'parametric'
To request the use of only the parametric solving scheme when solving a single first order ODE, possibly of high degree in dy/dx. Note however that by default dsolve tries to remove the parameter used during the solving process. To keep the parameter, specify the optional argument 'implicit' together with 'parametric'.

In the Examples section on that page you find an example (look for ode[4]);  you could also take a look at the link given there in the text.

Changing your input from vector to list makes your corrected code work:

restart;
z:=A*[x,y];
A:=1.0;
plots:-fieldplot(z,x=0..1,y=0..1);

Notice that all examples given in the help page for fieldplot has a list as the first argument. Notice also that the help page states:

A typical call to the fieldplot function is fieldplot(f, x=a..b, y=c..d), where f is a list of two expressions in x and y. a..b and c..d specifies the horizontal and vertical ranges of the field. The result is a two-dimensional vector field with the vector evaluated at (x, y) located at this point.
(emphasis added).

It is also mentioned that a VectorCalculus:-Vectorfield can be used, and indeed it works too:

restart;
z := VectorCalculus:-VectorField( A*<x, y>, 'cartesian'[x,y] );
A:=1.0;
plots:-fieldplot(z,x=0..1,y=0..1);
## Also this version works fine:

restart;
VectorCalculus:-BasisFormat(false);
z := VectorCalculus:-VectorField( A*<x, y>, 'cartesian'[x,y] );
A:=1.0;
plots:-fieldplot(z,x=0..1,y=0..1);

Try this:
First case:
 

restart;
ode:=diff(y(x),x,x)+diff(y(x),x)+y(x)=F;
bcs:=y(0)=1,y(2)=1;
F  := piecewise(x<=1, -x, x>1, 2*x+(x-1)^2);
sol:=dsolve({ode,bcs}); #That is not your solution
plot(rhs(sol),x=0..2); pE:=%:
plot(diff(rhs(sol),x),x=0..2);
plot(diff(rhs(sol),x,x),x=0..2,discont);
As:=dsolve({ode, bcs} ,numeric, maxmesh=500,abserr=1e-3); 
plots:-odeplot(As,[x,y(x)-rhs(sol)]); #Difference between numerical and exact solution: OK considering abserr

Now continue with the second case:

####
F := piecewise(x<=1, x^2+x+2, x>1, -x^2+x);
ode;
sol:=dsolve({ode,bcs});
plot(rhs(sol),x=0..2); 
plot(diff(rhs(sol),x),x=0..2);
plot(diff(rhs(sol),x,x),x=0..2,discont);
As:=dsolve({ode, bcs} ,numeric, maxmesh=4000,abserr=1e-3); 
plots:-odeplot(As,[x,y(x)-rhs(sol)]);

In the first case your exact "solution" on the interval 0..2 is less smooth than the solution provided by Maple's dsolve. Notice that if you rewrite the 2nd order equation as a first order system, then in your case one of the dependent variables (the one that was y'(x)) is not continuous. It should definitely be a requirement for a solution to a first order system on an interval a..b that the functions are (at least)  continuous on the whole interval a..b.
### To show the reason for this natural requirement, I shall show that if you give up continuity of the first derivative in your first case, then you have infinitely many solutions of your bvp problem:
 

F:= piecewise(x<=1, -x, x>1, 2*x+(x-1)^2);
sol1:=dsolve({ode,y(0)=1}) assuming x<1;
sol1:=subs(_C2=C1,sol1);
sol2:=dsolve({ode,y(2)=1}) assuming x>1;
sol2:=subs(_C2=C2,sol2);
eval(rhs(sol1),x=1)=eval(rhs(sol2),x=1); #One equation, two unknowns!
##So you can pick C2 to be anything, then C1 is determined.
##Solving for C1 in terms of C2:
solve(%,{C1});
sol11:=eval(sol1,%);
solT:=y(x)=piecewise(x<=1,rhs(sol11),rhs(sol2)); # "solution" on 0..2
plots:-animate(plot,[rhs(solT),x=0..2],C2=-1..1);
plots:-animate(plot,[diff(rhs(solT),x),x=0..2,discont],C2=-1..1);

The animation of y(x) above: