In a recent discussion about the existence of limit cycles I suggested using events using the mechanism you refer to:

http://www.mapleprimes.com/questions/207881-Convergence-Of-Nonlinear-Oscillator#answer223227

In reviewing an old worksheet of mine I found the following example (competing species). That should give the general idea.

restart;
f:=(x,y)->x*(e1-s1*x-a1*y);
g:=(x,y)->y*(e2-s2*y-a2*x);
#The system:
sys:={diff(x(t),t)=f(x(t),y(t)),diff(y(t),t)=g(x(t),y(t))};
#The system sys has 6 parameters
#We non-dimensionalize:
varchange:={x(t)=X*xi(tau),y(t)=Y*eta(tau),t=T*tau}; #X,Y,T constants to be chosen
PDEtools:-dchange(varchange,sys,[tau,xi,eta]);
solve(%,{diff(xi(tau),tau),diff(eta(tau),tau)});
systemp:=expand(%);
#Now we have lots of choices possible.
#One is to make the coefficient of xi in the ode for xi equal to 1:
eval(systemp,{T=1/e1});
#By choosing X=e1/s1 and Y=e1/a1 we get rid of all parameters in the ode for xi: just a choice.
#So doing all 3 at once:
eval(systemp,{T=1/e1,X=e1/s1,Y=e1/a1});
#We may give the coefficents in the ode for eta new names:
sys_nondim:=subs(a2=a*s1,s2=s*a1,e2=e*e1,%);
##
##So we have the system
    {diff(eta(tau), tau) = -xi(tau)*eta(tau)*a-eta(tau)^2*s+eta(tau)*e, diff(xi(tau), tau) = -xi(tau)^2-        xi(tau)*eta(tau)+xi(tau)}

It follows from the ode for xi that xi and xi^2 have the same dimension. Thus xi is dimensionless. By the same ode then also eta is dimensionless and further that tau must be.
From the ode for eta it then follows that the new coefficients: a, s, and e are also dimensionless. These new coefficients (a=a2/s1, s=s2/a1, e=e2/e1) are sometimes referred to as dimensionless groups (or some such term).
It is noteworthy that the number of parameters in sys_nondim is only 3 (instead of the original 6). This makes it easier to analyze the system.

N:=5:
K:=Matrix(N+1,(i,j)->D[1$(i-1),2$(j-1)](k)(0,0)/(i-1)!/(j-1)!);

or you could make K a function of N:

restart;
K:=N->Matrix(N+1,(i,j)->D[1$(i-1),2$(j-1)](k)(0,0)/(i-1)!/(j-1)!);
##Testing:
K(3);

Your second initial condition should use D:
bw1:=D(r)(0) = 0;

res:=convert(taylor(y(x+h), h=0, 3), polynom);
params:={h=.1,y(x)=7,D(y)(x)=9,(D@@2)(y)(x)=13};
eval(res,params);

I copied you line:

font_size := thickness = 4, font = [bold, "TimesNewRoman", 30], labelfont = [bold, "TimesNewRoman", 30], axis[1] = [thickness = 3], axis[2] = [thickness = 3], captionfont = [bold, "TimesNewRoman", 30], legendstyle = [font = [bold, "TimesNewRoman", 30]];

into a worksheet in Maple 2015.2. I exported it as an mpl file using the menu.

Then I tried two different approaches:
1. I opened the file in Maple just as I usually open worksheets, but choosing file type as mpl in the "Open" window.
This worked without any problem.
2. I tried reading it in using the read statement, in my case:
  read "E:/MapleDiverse/MaplePrimes/testFONT.mpl";
That worked too.

So I wonder what you did.

You should not have multiplication signs after tanh and sech. So it should be

u := -3*sqrt(mu)*tanh(A+sqrt(-mu)*(x-(1/2)*mu*t))/sqrt(6*a)+3*sqrt(mu)*sech(A+sqrt(-mu)*(x-(1/2)*mu*t))/sqrt(-6*a);

You have a difference of two huge positive terms when x is just moderately large.

restart;
S1:=Sum(pochhammer(1/3,k)*3^k*x^(3*k)/(3*k)! ,k=0..infinity);
S2:=Sum(pochhammer(2/3,k)*3^k*x^(3*k+2)/(3*k+2)!  ,k=0..infinity);
S3:=Sum(pochhammer(2/3,k)*3^k*x^(3*k+1)/(3*k+1)!  ,k=0..infinity);
S4:=Sum(pochhammer(1/3,k)*3^k*x^(3*k+1)/(3*k+1)!  ,k=0..infinity);
s1,s2,s3,s4:=op(value~([S1,S2,S3,S4]));
## So you are considering:
s1*s2 - s3*s4;

#But just try plotting these 4 (and notice I only let x go to 10):
plot([s1,s2,s3,s4],x=0..10,legend=["s1","s2","s3","s4"],thickness=[1,1,3,3]);
## and the inert ones, but set Digits high:
Digits:=50:
plot([S1,S2,S3,S4],x=0..10);
## You notice from the first plot that it seems not unlikely that s1*s2 is of the same order of magnitude as s3*s4. Thus you have a potentially severe numerical problem because these quantities become huge.
####Note about asympt: This won't work unless additional facilities are added. So skip to the last 5 lines:
You could look at the asymptotic behavior of these:
asympt~([s1,s2,s3,s4],x,6);
A:=convert(%,polynom);
plot(A[1]*A[2]-A[3]*A[4],x=10..18);
#This seems to confirm the statement about huge cancellations in subtraction.
####
If you could rewrite the difference S1*S2-S3*S4 to avoid this problem, you wouldn't have to raise Digits to enormous values.
Try this:
evalf[50](eval(s1*s2-s3*s4,x=100));
evalf[500](eval(s1*s2-s3*s4,x=100));
evalf[1000](eval(s1*s2-s3*s4,x=100));

You should use unapply in f:

restart;
with(CodeGeneration):
f := proc(r) unapply(r,x) end proc;
to_translate := f(convert(series(sin(x),x,20),polynom));
##Test
to_translate(7);
CodeGeneration['C']( to_translate );

## Compare with your version:
restart;
f := proc(r) x->r end proc;
to_translate := f(convert(series(sin(x),x,20),polynom));
## Test
to_translate(7);

## Another version that works:
restart;
f := proc(r) subs(rr=r,x->rr) end proc;
to_translate := f(convert(series(sin(x),x,20),polynom));
to_translate(7);

Here is a way:

psi := x^(4/5)*f(x, eta)/(1+x)^(1/20);
psi1:=eval(psi,eta = y/(x^(1/5)*(1+x)^(1/20)));
diff(psi1, x);
diff(psi1, y);

restart;
with(inttrans);
eq1:=(1/2)*x^2+(1/2)*x^3+(1/12)*x^4 = int((x-t-1)*u(t)+(x-t+1)*v(t), t = 0 .. x);
IntegrationTools:-Expand(eq1);
eq1s:=laplace(%,x,s);
eq2:=(3/2)*x^2-(1/6)*x^3+(1/12)*x^4 = int((x-t+1)*u(t)+(x-t-1)*v(t), t = 0 .. x);
IntegrationTools:-Expand(eq2);
eq2s:=laplace(%,x,s);
##
solve({eq1s,eq2s},{laplace(u(x),x,s),laplace(v(x),x,s)});
invlaplace(%,s,x);

  

Don't ever try to execute a whole block of statements without having executed them one after the other and observed the output.
I split your whole block of code into its single commands and started executing from the top. I didn't get very far. After having executed Digits:=13, executing the next one:

sys_ode := diff(ph(t), t) = urj*prj(t)+urd*prd(t)+ugd*pgd(t)-yc*ph(t), diff(pc(t), t) = yc*ph(t)-(yrj+yrd+ygd+yf)*pc(t), diff(prj(t), t) = yrj*pc(t)-urj*prj(t), diff(prd(t), t) = -urd*prd(t)+yrd*pc(t), diff(pgd(t), t) = -ugd*pgd(t)+ygd*pc(t), diff(pf(t), t) = yf*pc(t)*simplify(add(u, u = sys_ode));

resulted in the error message
   
Error, recursive assignment

That error is issued since you have the as yet unassigned sys_ode on the left of := but also on the right buried in add(u, u = sys_ode).

This assignment would make sense, if sys_ode had previously been defined to be something, but it hasn't.
A very simple example of the same thing:

restart;
a:=a+1;
Error, recursive assignment

Do you mean that by R = 2*b*(b^2-9*b-6)/(10*b^3+21*b^2+16*b+4) is defined a function of one variable (b)?

If so, you will notice from the plot below (and in other ways) that R is negative for b=1..5 and not particularly large in absolute value.

restart;
R := 2*b*(b^2-9*b-6)/(10*b^3+21*b^2+16*b+4);
plot(R,b=1..5);
solve(R=0,b);
evalf(%);

restart;
E:=-2*P1*A[Q1]+A[0]^2+A[Q1]^2-2*A[0]*rho*P2/(rho+Q1)+A[Q3]^2-2*P3*A[Q3]+P1^2+rho^2*P2^2/(rho+Q1)^2+P3^2+m^2*c^4;
mtaylor(E,[A[Q1]=P1,A[0]=P2*rho/(rho+Q1),A[Q3]=P3],6);
## You can obtain the same result with:
Student:-Precalculus:-CompleteSquare(E,[A[0],A[Q1],A[Q3]]);

In your second example (a first order system of 8 odes) you have a function, which basically switches from omega[sw] to omega[st] in a very short time. The way you have implemented it, the switch is continuous. I'm assuming that this switch is implemented in an ad hoc manner and could as well be done differently.

Thus I shall present a way to do it discontinuosly by using the events option in dsolve/numeric:

restart;
Digits:=15:
K:=Matrix([<0, -1, 1, -1>,<-1, 0, -1, 1>,<-1, 1, 0,-1>,<1, -1, -1,0>]);
for i to 4 do
  r:=sqrt((u[i](t))^2+(v[i](t))^2);
  Equ[i]:=diff(u[i](t),t)=Au*(1-r^2)*u[i](t)-omega[i](t)*v[i](t);
  Eqv[i]:=diff(v[i](t),t)=Av*(1-r^2)*v[i](t)+omega[i](t)*u[i](t)+(K. <seq(v[i](t),i=1..4)>)[i];
  EqSys[i]:=[Equ[i],Eqv[i]]
end do:
params:={omega[st]=2*Pi,omega[sw]=4*2*Pi,Au=5,Av=5};

sys:={seq(EqSys[i][],i=1..4)};
##In ic the boolean variable b[i](t) is given initial value 0 if v[i](0) < 0 and 1 otherwise:
ic:={u[1](0)=0, v[1](0)=-.01,u[2](0)=0, v[2](0)=-0.1,u[3](0)=0, v[3](0)=0.1,u[4](0)=0, v[4](0)=0.1,b[1](0)=0,b[2](0)=0,b[3](0)=1,b[4](0)=1};
## The b[i]'s take on values 0 and 1 only:
sysE:=subs(seq(omega[i](t)=b[i](t)*omega[st]+(1-b[i](t))*omega[sw],i=1..4),sys);
resE:=dsolve(eval(sysE union ic,params), numeric,maxfun=0,discrete_variables=[seq(b[i](t)::boolean,i=1..4)],
        events=[seq([side(v[i](t),post),b[i](t)=1-b[i](t)],i=1..4)]);
plots:-odeplot(resE,[u[1](t),v[1](t)],0..15,numpoints=2000);
plots:-odeplot(resE,[u[1](t),v[1](t)],100-.6..100,numpoints=2000);
plots:-odeplot(resE,[[t,u[1](t)],[t,v[1](t)]],0..15,numpoints=2000,size=[1800,200]);
plots:-odeplot(resE,[t,u[1](t)],0..15,numpoints=2000,size=[1800,200]);
plots:-odeplot(resE,[[t,b[1](t)],[t,v[1](t)]],0..15,numpoints=2000,size=[1800,200],thickness=[3,1]);
plots:-odeplot(resE,[u[1](t),v[1](t),u[2](t)],0..15,numpoints=2000);
plots:-odeplot(resE,[u[1](t),v[1](t),u[3](t)],10..15,numpoints=2000);
##For an explanation of the trigger side(y(t),post), see:
?dsolve,numeric,events