I think you need to contact Maplesoft regarding the licensing. After that do a fresh installation on the new machine. Don't just copy from a folder on the old machine to the new.

The animation parameter in animate will be made a floating point number during execution, thus an expression intended to be sols[1..5] will be seen as sols[1..5.0], which won't work. The remedy is simple: use round as in this example, where I have used the example points from Kitonum's answer:

sols := [[0,0], [6,0], [6,4], [0,4], [3,6], [6,4], [0,0], [0,4], [6,0]];
plots:-animate(plot,[sols[1..round(i)],thickness=3],i=1..nops(sols),paraminfo=false);

The code you sent me was:
 

restart; 
with(plots): with(VectorCalculus): with(LinearAlgebra):
addcoords(EEspherical, [rho, theta, phi], [rho*cos(phi)*sin(theta), rho*sin(phi)*sin(theta), rho*cos(theta)]);

P := proc (rho, theta, phi) options operator, arrow; `<,>`(rho*cos(phi)*sin(theta), rho*sin(phi)*sin(theta), rho*cos(theta)) end proc;
 
`vel&rho;` := (D[1](P))(rho, theta, phi); 
`vel&theta;` := (D[2](P))(rho, theta, phi); 
`vel&phi;` := (D[3](P))(rho, theta, phi);
 
Norm(`vel&theta;`, 2);
## Now just continue with
simplify(%) assuming real;

Alternatively, you could do:

Norm(`velθ`, 2,conjugate=false);
simplify(%) assuming rho>0;

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);