Besides all the occurrences of [ ] where ( ) should be used (quite a few) as Carl pointed out, there is one place where you use { } instead of ( ). The curly braces { } are only used for sets. So correct that in the section Taxas, where r[1] is defined.
If you add maxfun=10^5 to the dsolve command plotting is no problem.
Here is the graph of T:

Using the form introduced by vv in his answer.

 

restart;
J:=Int( ln(1-y)*y^a*(1-y)^b, y=0..1);
Jb:=applyop(int,1,J,b=0..b); #Physicists' abuse of notation OK in recent versions
Jb:=eval(applyop(int,1,J,b=0..bb),bb=b); # My preference
res1:=value(Jb) assuming b>-1,a>-1;
res2:=diff(res1,b); # Result

res2 := Psi(b+1)*GAMMA(b+1)*GAMMA(a+1)/GAMMA(a+b+2)-GAMMA(b+1)*GAMMA(a+1)*Psi(a+b+2)/GAMMA(a+b+2)

## Actually, we might as well use indefinite integration:

Jbindef:=applyop(int,1,J,b);
res1a:=value(Jbindef) assuming b>-1,a>-1;
res3:=diff(res1a,b);

The result is in another form:

res3 := (Psi(b+1)-Psi(a+b+2))*Beta(a+1, b+1)

But
simplify(convert(res2-res3,GAMMA)); 
returns 0.

I had success with the first four of the HAs using continuation:

for k to nops(HA) do
    res[k] := dsolve(eval({Eq1, Eq2, Eq3, Eq4, IC1, IC2}, params union {K1 = lambda*HA[k]}), numeric,maxmesh=1024,continuation=lambda)
end do;

But not immediately for the fifth, i.e. the value 0.8.

It appears that InitialValueProblem wants diff(y(t),t) instead of y'(t) or D(y)(t) (as the prime version is interpreted in 2D-math input).
So use:
InitialValueProblem(diff(y(t),t) = t*y(t)+1/y(t), y(0) = 3, t = 2, method = euler, numsteps = 5, output = solution);

I modified the example you gave because the function doesn't have any critical points.
 

restart;
f:=-2*x^2+x+3*y-y^3;
plot3d(f,x=-5..5,y=-5..5);
## Finding the critical points:
eqs:={diff(f,x)=0,diff(f,y)=0};
res:=solve(eqs,{x,y});
pts:=map(subs,[res],[x,y]);
## Simplest is it to use this package:
with(Student:-MultivariateCalculus):
SecondDerivativeTest(f,[x,y]=pts);
SecondDerivativeTest(f,[x,y]=pts,output=hessian);
## More on your own you can do this instead:
H:=VectorCalculus:-Hessian(f,[x,y]);
H1:=eval(H,res[1]);
LinearAlgebra:-Determinant(H1);
LinearAlgebra:-Trace(H1);
H2:=eval(H,res[2]);
LinearAlgebra:-Determinant(H2);

Since nobody has mentioned fnormal in this context, I will call attention to it:

restart;
fnormal(evalf(Pi-3.14),10,1e-2);
fnormal(evalf(Pi-3.14),10,1e-3);
fnormal(evalf(Pi-3.14),4);
fnormal(evalf(Pi-3.14),5);
# I have kernelopts(floatPi=false) in my maple.ini file.

Correct the syntax as I mentioned, Use a high setting of Digits, perhaps 50.
Change the lower end of range to ep in dsolve and odeplot.
Digits:=50:
ics:={ P(ep) = 1.3668*10^14, m(ep) = 0, rho(ep)=4.2983*10^(14),v(ep)=-0.7342};
sol := dsolve(sys union ics, numeric,method=rosenbrock, range = ep .. 6*10^5); # Upper end changed since you won't get to 20*10^5
## You will have a problem at about r = 567836.24.
## Basically because rho(r) will be very near to zero and the numeric solver will have problems with e.g. rho(r)^(2/3)

plots:-odeplot(sol,[r,rho(r)],ep..567837,thickness=3);

plots:-odeplot(sol,[r,rho(r)],500000..567837,thickness=3);
############ A simple example that shows the problem:

restart;
ode:=diff(rho(r),r)=-rho(r)^(2/3);
res:=dsolve({ode,rho(0)=1},numeric);
res(5); #Gives the error:

Error, (in res) cannot evaluate the solution further right of 2.9995012, probably a singularity

## Plotting

plots:-odeplot(res,[r,rho(r)],0..5);
You get a warning:

Warning, cannot evaluate the solution further right of 2.9995012, probably a singularity
## This is not just some quirk in the numerics:

sol:=dsolve({ode,rho(0)=1});
plot(rhs(sol),r=0..5);
## However:
rest:=odetest(sol,ode);
simplify(rest) assuming r<3;
simplify(rest) assuming r>3;
## A further comment: At r=3 the solution can be extended as identically zero for r>3.

########## Possible solution here in this simple example:

restart;
ode:=diff(rho(r),r)=-piecewise(rho(r)>0,rho(r)^(2/3),0);
res:=dsolve({ode,rho(0)=1},numeric);
res(5);
plots:-odeplot(res,[r,rho(r)],0..5);
#####################
The idea from the simple example works.
I bring all the code:
 

restart;
G := 6.6743*10^(-8); 
a := 2.1197*10^24; 
b := 1/3.035; 
c := 2.99792458*10^10; 
d := 2.035; 
pi = 3.143; # Got nothing to do with Pi, I assume?
polytrop := 5/3; 
K := 5.38*10^9/(c^2*(G/c^2)^(2/3));
sys := {diff(P(r), r) = -G*(a*P(r)^b+(1+d)*P(r)/d)*(m(r)+4*Pi*r^3*P(r)/c^2)/(c^2*(r^2-2*G*r*m(r)/c^2)), diff(m(r), r) = 4*Pi*rho(r)*(1+K*rho(r)^(polytrop-1)/(polytrop-1))*r^2, diff(rho(r), r) = -(rho(r)*(1+K*rho(r)^(polytrop-1)/(polytrop-1))+K*rho(r)^polytrop)*(4*Pi*r^3*K*rho(r)^polytrop+m(r))/(K*polytrop*rho(r)^(polytrop-1)*r*(r-2*m(r))), diff(v(r), r) = 1.485232054*10^(-28)*(m(r)+4.450600224*10^(-21)*Pi*r^3*P(r))/(r^2-1.485232054*10^(-28)*r*m(r))};
RHO:=select(has,sys,diff(rho(r),r));
M:=select(has,sys,diff(m(r),r));
rho_eq:=diff(rho(r),r)=piecewise(rho(r)>0,rhs(op(RHO)),0);
m_eq:=diff(m(r),r)=piecewise(rho(r)>0,rhs(op(M)),0);
newsys:=sys minus (RHO union M) union {rho_eq,m_eq};
ep := 0.1e-10;
ics:={ P(ep) = 1.3668*10^14, m(ep) = 0, rho(ep)=4.2983*10^(14),v(ep)=-0.7342};
Digits:=50:
sol := dsolve(newsys union ics, numeric,method=rosenbrock, range = ep .. 2*10^6);
plots:-odeplot(sol,[r,rho(r)],ep..2*10^6,thickness=3);
plots:-odeplot(sol,[r,m(r)],ep..2*10^6);
plots:-odeplot(sol,[r,P(r)],ep..2*10^6);
plots:-odeplot(sol,[r,v(r)],ep..2*10^6);

Example plot. The plot of m:

You have u*(x, 0) = 0 .
It should be u(x,0) = 0. If you are using 2D input then remember that a space is interpreted as *.
If you think that is horrible, then switch to 1D math input, aka Maple Input. Go to Yools/Options/Display/Input Display.

Often there may be a parameter involved, so animation may be relevant.
Try the following code several times without restart. The matrix will change and so will the animation.
 

with(LinearAlgebra):
A:=Matrix(4,(i,j)->if i=2 and j=3 then a else rand(-9..9)() end if);
L:=Eigenvalues(A):
plots:-animate(plots:-complexplot,[convert(L,list),style=point,symbolsize=20,symbol=solidcircle],a=-5..5);

I tried the following:

stopat(s1,4);
Then I clicked on the button next.
After that I entered
Y:=eval(y);
in the space provided for debugger commands.
Then pressed the button Execute followed by pressing quit.
In the worksheet I then did
eval(Y);
You may want to finish with
unstopat(s1);
 

fsolve only returns one solution for an expression that is not a polynomial.
See ?fsolve
But there are other ways:

fsolve(sin(x),x=1..17);
RootFinding:-Analytic(sin(x),x=1-.1*I..17+.1*I);
sort([%]);
identify~(%);
solve({sin(x),x>1,x<17},allsolutions,explicit);

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;