I checked with one much older version Maple 12 as well as with Maple 2021.1.

As far as the display of the leading minus is concerned the behavior is the same.
Observe the output from this code:
 

restart;
t__FET_TurnOff := solve(-C__iss2*V__gs_th2 + Q__total = i__GateDriveN*t__off, t__off);

numer(t__FET_TurnOff);
denom(t__FET_TurnOff);
%%/%;
(numer/denom)(t__FET_TurnOff); 

To show numerical evidence of convergence of solutions to a single (periodic) solution try this:
 

restart;
ode:=diff(x(t),t)=-(1+1/10*sin(10*t))*x(t)+1;
res:=dsolve({ode,x(t0)=x0},numeric,parameters=[t0,x0]);
res(parameters=[0,1]); #setting parameters
plots:-odeplot(res,[t,x(t)],t=0..20,size=[1200,default],numpoints=5000); p1:=%:
res(parameters=[0,2]);
plots:-odeplot(res,[t,x(t)],t=0..20,size=[1200,default],numpoints=5000,color=blue,view=0.9..1.1); p2:=%:
plots:-display(p1,p2);
res(parameters=[0,-1]);
plots:-odeplot(res,[t,x(t)],t=0..20,size=[1200,default],numpoints=5000,color=cyan,view=0.9..1.1); p3:=%:
plots:-display(p1,p2,p3);
res(parameters=[2,-1]);
plots:-odeplot(res,[t,x(t)],t=0..20,size=[1200,default],numpoints=5000,color="Green",view=0.9..1.1); p4:=%:
plots:-display(p1,p2,p3,p4);
plots:-display(p1,plot(sin(10*t+1.7)/100+1,t=0..20,color=blue));
Q:=proc(t0,x0,{range::range:=0..10,numpoints::posint:=5000,size::[posint,posint]:=[1200,default]}) if not [t0,x0]::list(realcons) then return 'procname'(_passed) end if;
   res(parameters=[t0,x0]);
   plots:-odeplot(res,[t,x(t)],range,:-numpoints=numpoints,:-size=size,_rest)
end proc;
Q(0,2,range=0..20,view=0.9..1.1);
plots:-animate(Q,[0,x0],x0=-2..3,trace=24,view=0.9..1.1,size=[1200,400]);

The animation at the end:

Your ode is very similar to the one in https://mapleprimes.com/questions/232633-How-To-Solve-Nonlinear-Ode-With-Numeric-Method

to which I responded. Your ode can similarly be integrated once to result in
 

ode := delta*diff(u(y), y)^3 + diff(u(y), y) = p*y

where use of D(u)(0) = 0 has already been done. Proceeding as in the link in this case leads first to 3 odes, but 2 of these can be eliminated since they don't satisfy D(u)(0) = 0 ( assuming delta>0).
The remaining is

diff(u(y), y) = 1/6*12^(1/3)*((9*p*y*sqrt(delta) + sqrt(27*delta*p^2*y^2 + 4)*sqrt(3))^(2/3) - 12^(1/3))/(sqrt(delta)*(9*p*y*sqrt(delta) + sqrt(27*delta*p^2*y^2 + 4)*sqrt(3))^(1/3))

Which leads to the solution by direct integration using that u(k) =1:
 

u(y) = Int(1/6*(2*18^(1/3)*((9*p*_z1*sqrt(delta) + sqrt(81*_z1^2*delta*p^2 + 12))^2)^(1/3) - 12)/(sqrt(delta)*(108*p*_z1*sqrt(delta) + 12*sqrt(81*_z1^2*delta*p^2 + 12))^(1/3)), _z1 = k .. y) + 1

Unfortunately this integral cannot be done symbolically unlike the case considered in the link.

I know you want ode solved numerically and also by a shooting method.

Here though, I shall point out that your ode bvp can be solved exactly.

## Your worksheet has a syntactical error in bc:  You have D*u(0) where it should be D(u)(0).

restart;
#p:=2:
ode := diff(u(y), y, y) + beta*diff(diff(u(y), y)^2, y) = p;
bcs:=D(u)(0)=0,u(h)=1;
odeC:=map(int,ode,y)+(0=C);
eval[recurse](convert(odeC,D),{bcs[1],y=0});
ode1:=eval(odeC,C=0);
solve(ode1,{diff(u(y),y)});
odes:=op~([%]);
eval[recurse](convert(odes,D),{y=0,bcs[1]});
# So only odes[1] satisfies bcs[1]
sol1:=dsolve({odes[1],u(h)=1});
# Dealing with the original second order ode:
sol:=dsolve({ode,bcs});
simplify(sol-sol1); # OK

For a numerical solution you also need to provide a numerial value for h.

I don't see how you can do shooting on the first order ode (odes[1])  which already satisfies D(u)(0)=0.
But obviously you can use a numerical method as in the following:
 

p:=2;
res:=dsolve({odes[1],u(h)=1},numeric,parameters=[beta,h],abserr=1e-12,relerr=1e-9);
res(parameters=[beta=1/8,h=3]);
plots:-odeplot(res,[y,u(y)],0..3);
plots:-odeplot(res,[y,u(y)-eval(rhs(sol1),{beta=1/8,h=3})],0..3);#plot of error

Very good agreement between the exact and numerical solutions for h=3, beta=1/8.

######### In a reply below I have added a link to an expanded version of the code above.

It includes conversion of the given ode to a traditional system of 2 odes and has also been solved by shooting for that system.

For the 2 matrices (or more) you can just use <... , ... > in this way:
 

A:=LinearAlgebra:-RandomMatrix(3,2);
B:=LinearAlgebra:-RandomMatrix(4,2);
C:=LinearAlgebra:-RandomMatrix(2,2);
M:=<A,B>;
M1:=<A,B,C>;
M2:=<M,C>;
M1-M2;

This works:
 

restart;

expr:=x*A[2]+A[1];
select(type,expr,identical(x) &* anything);
select(type,expr,identical(x^2) &* anything);

See ? type, structured

Using simple shooting shows that there can be several solutions.

I take here the example where k = 1.5, alpha = 2, and infinity = 10.
In this case we find 3 solutions, 2 of which may be considered unphysical or maybe rather not corresponding properly to an acceptance of infinity = 10.

restart;
odeSys := f(xi)+3*diff(g(xi), xi)-xi*diff(f(xi),xi)=0,
            f(xi)^2+3*g(xi)*diff(f(xi),xi)-xi*f(xi)*diff(f(xi),xi)=3*diff(f(xi),xi,xi)+18*k*diff(f(xi),xi)^2*diff(f(xi),xi,xi);
  bcs := f(0)=alpha, f(infinity)=0, g(0)=0;
inf:=10;
respar:=dsolve(eval({odeSys,f(0)=alpha,g(0)=0,D(f)(0)=f1},{alpha=2,k=1.5}),numeric,parameters=[f1],output=listprocedure);
F:=subs(respar,f(xi));
Q:=proc(f1) if not f1::realcons then return 'procname(_passed)' end if;
   respar(parameters=[f1]);
   F(inf)
end proc;
plot(Q(f1),f1=-2..0);
f1val1:=fsolve(Q,-1.8..-1.6);
f1val2:=fsolve(Q,-1..-0.98);
f1val3:=fsolve(Q,-0.95..-0.9);
f1vals:=[f1val1,f1val2,f1val3];
for i from 1 to numelems(f1vals) do
  Q(f1vals[i]);
  p[i]:=plots:-odeplot(respar,[xi,f(xi)],0..inf,caption=typeset(D(f)(0)=f1vals[i]))
end do:
plots:-display(seq(p[i],i=1..numelems(f1vals)),insequence);
plots:-display(seq(p[i],i=1..numelems(f1vals)),caption=typeset(D(f)(0)=f1vals));

Notice that Q sets the parameter f1 each time it is called.
The the 3 solutions in an animation:

Comment.
Making the procedure Q have 2 arguments, f1 and inf (now unassigned, was 10) and trying much higher values of inf (e.g. 100) shows that 3 solutions of Q = 0 remain. They are all negative. Also it appears that the smallest in absolute value is roughly the same for all values of inf. The middle one moves ever closer to the first while the absolute value of the third root becomes rather large.
Values for inf=100:  -363.8130077, -0.9332179835, -0.9272464476.

My guess is that as inf -> infinity the 3rd root tends to -infinity and the first and second root becomes a double root. Thus we are left with one root.

I tried your worksheet in Maple 18, Maple 2020.2, and in Maple 2021.1.

The double loop over i and j got stuck in Maple 2020.2, but not in the other two.

Already the computation of evalf(int(r*h0[i]*g[j], r = 0 .. 1)) gets stuck for i = j = 1.
Try this version instead:
 

for i to n do
    for j to n do h[j][i] := int(r*h0[i]*g[j], r = 0 .. 1,numeric,method=_d01ajc)/int(r*(-r^2 + 1)*g[j]^2, r = 0 .. 1,numeric,method=_d01ajc); end do;
end do;

The use of method=_d01ajc makes the computation real fast.

Try this:
 

expr:=piecewise(x(t)<0, -1, x(t)<1, 0, 1) + f(t) + piecewise(t<0, 0, t<1, 1, 0) + x(t)*piecewise(t<0, 1, t<1, 0, 1) + a*piecewise(t<0, 3, t<1, -1, 2);
####
select(hastype,expr,specfunc(piecewise));
select(hastype,%,And(relation,satisfies(r->lhs(r)=t)));

Trying this:

plot(w->int(eval(A1+A2,omega=w),theta=0..2*Pi,numeric,digits=5),0..50);

I got

You won't be able to find a symbolic solution. You must solve numerically.

Thus you must give values to the parameters c and t.
You must also give boundary or intial values.

I'm assuming that by e^() you mean the exponential function, which in Maple is exp. So I have made that replacement.
Here is an example of solving the initial value problem.

restart;
sys := -diff(g(x), x)*exp(c*x)*c/(3*t) + f(x) = 0, t*(c*exp(2*c*x)*diff(g(x), x) + t*exp(c*x)*(2*c^2 + 3*f(x) - 1))*diff(f(x), x) + 6*c^2*(diff(f(x), x, x)*t*exp(2*c*x) + diff(f(x), x, x, x)*exp(3*c*x)/3) = 0;
sol:=dsolve({sys}); # essentially unsolved, but work has been done
dsolve(sol[1]); # returns NULL i.e. nothing
#Solving the original system numerically:
res:=dsolve({sys,f(0)=f0,D(f)(0)=f1,(D@@2)(f)(0)=f2,g(0)=g0},numeric,parameters=[c,t,f0,f1,f2,g0]);
#Setting parameters (an example):
res(parameters=[1,2,0.5,0.6,0,0]);
plots:-odeplot(res,[[x,f(x)],[x,g(x)]],-1..10);
## Here is a procedure that can be used by animate:
Q:=proc(c,t,f0,f1,f2,g0,scene::list:=[[x,f(x)],[x,g(x)]],x_interval::range:=-1..10) 
   if not [c,t,f0,f1,f2,g0]::list(realcons) then return 'procname(_passed)' end if;
   res(parameters=[c,t,f0,f1,f2,g0]);
   plots:-odeplot(res,scene,x_interval,_rest)
end proc;
## Examples of use all with c, t, f1, f2, g0 fixed
plots:-animate(Q,[1,2,f0,0.6,0.1,-.3,0..20],f0=-1..1,trace=24); 
plots:-animate(Q,[1,2,f0,0.6,0.1,-0.3,[f(x),diff(f(x),x)],0..30],f0=-1..1);
plots:-animate(Q,[1,2,f0,0.6,0.1,-0.3,[seq([x,diff(f(x),[x$k])],k=0..2)],-0.3..5],f0=0..5);

In the worksheet attached I end up displaying all the separatrices in sequence producing the animation below.
MaplePrimes21-05-05_separatrices.mw

Since the problem is formulated in vector form the Maple solution should be given as such, I think.
Even Maple 2021 doesn't quite do that.
Here is the closest I got:
 

restart;
X:=<x1(t),x2(t)>;
Sys:=diff(X,t)=<-4,-2;6,3>.X+<2/(exp(t)+1),-3/(exp(t)+1)>; 
ics:=eval(X,t=0) = <0,ln(2)>;
 
sol1:= dsolve([ Sys, ics]); # Not on vector form
sol:=X=eval(X,sol1); # Now on vector form

The solution looks like this:

You have a classical Volterra-Lotka system.
All solutions are periodic. The orbits in phase space are closed.
Use odeplot with a numerical solution as done here.
 

restart;
sys := diff(x1(t), t) = 3*x1(t) - 2/1000*x1(t)*x2(t), diff(x2(t), t) = 6/10000*x1(t)*x2(t) - 1/2*x2(t);
ics := x1(0) = 1000, x2(0) = 500;
res:=dsolve({sys,ics},numeric);
plots:-odeplot(res,[[t,x1(t)],[t,x2(t)]],0..15,legend=[x1,x2],size=[1000,500]);
plots:-odeplot(res,[x1(t),x2(t)],0..7,frames=100);

The last is an animation of the orbit in phse space and is seen below.

I don't have Maple 13, but I have Maple 12. It gives the same error.
The error is NOT related to the way you produce the subscript, which is seen by this code:
 

restart;
H1:=v*P[r]-((1-P[r]^2*r/(P[r]^2*r-1))/r)^(1/2);
simplify(H1); ## Error in Maple 12
normal(H1);   ## OK
radnormal(H1);         ## Same error
simplify(H1,symbolic); ## Same error

The error is also present in Maple 15, 16, 17, 18, and 2015, but not in Maple 2016.