Here is the code for using dsolve/numeric and odeplot. No need to write the ode as a first order system. dsolve will actually do that itself.
 

ode:=diff(u(t),t,t)-mu*(c-b*u(t)^2)*diff(u(t),t)+u(t)=f*cos(omega*t);
params:={mu=0.001,b=2.4,c=2.4,f=8.28,omega=0.182};
res:=dsolve({eval(ode,params),u(0)=3,D(u)(0)=4},numeric,range=0..500,maxfun=10^5);
plots:-odeplot(res,[u(t),diff(u(t),t)],100..500,refine=1);

Notice that by using the interval 100..500 in odeplot we won't be seeing any transient behavior.
The use of the range argument in dsolve allows us to use refine in odeplot. The alternative is to use a high value of numpoints (like 10000).
By keeping parameters as equations in a variable I called params we can easily change those parameters.
I get this picture, which doesn't reproduce the picture in the reference:

Here follows an animation in mu. I have removed mu from params and used the parameters option in dsolve.
 

params1:={b=2.4,c=2.4,f=8.28,omega=0.182};
resP:=dsolve({eval(ode,params1),u(0)=3,D(u)(0)=4},numeric,range=0..500,maxfun=10^5,parameters=[mu]);
pp:=proc(mu) resP(parameters=[mu]); plots:-odeplot(resP,[u(t),diff(u(t),t)],100..500,refine=1) end proc;
plots:-animate(pp,[mu],mu=0.001..0.1,frames=50);

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?