The two problems give different answers, but that is not surprising, since you have two quite different odes.

The first one is in fact:
fx:=int(exp(-(1/10)*t^2), t = 0 .. x);
dsolve({diff(y(x), x) = y(x)+fx, y(0) = 0}); # Symbolic solution for comparison with your numeric:
evalf(eval(%,x=1));
                    
The second is in fact:
f2:=eval(fx,x=1);
dsolve({diff(y(x), x) = y(x)+f2, y(0) = 0}); # Symbolic solution for comparison with your numeric:
evalf(eval(%,x=1));
                     
These results agree with your numeric results.

Only you can decide which of the two differential equations is relevant for you to solve.

Replace AA by:
AA2:=subs(Y=Y(r1,r2,r3),AA);

so the dependency on r1,r2,r3 becomes explicit. Then dchange works.

restart;
for i to 4 do
  v:=L/2*cos(w*t+phi[i]);
  Eq[i]:=u[i](t)=v
end do;

or even simpler

restart;
for i to 4 do   Eq[i]:=u[i](t)=L/2*cos(w*t+phi[i]) end do;

You can use select:

Pz:=select(s->vdot(s)>0,zeros);
##If you want the nonnegative ones too, you can do both at the same time:
Pz,Nz:=selectremove(s->vdot(s)>0,zeros);

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