I tried the following using worksheet interface and Maple input ( as always):

restart;
U:=-(s^2/(s^2+1)+1/(s^2+1)-1)*exp(sqrt(s)*x/sqrt(phi))/(s^2+1)+exp(-sqrt(s)*x/sqrt(phi))/(s^2+1):
inttrans:-invlaplace(U,s,t) assuming phi>0,x>0;
                       
#Whereas the following, where U(x) is a remember table assignment to the function U:
restart;
U(x):=-(s^2/(s^2+1)+1/(s^2+1)-1)*exp(sqrt(s)*x/sqrt(phi))/(s^2+1)+exp(-sqrt(s)*x/sqrt(phi))/(s^2+1):
inttrans:-invlaplace(U(x),s,t) assuming phi>0,x>0;
                         U(x) Dirac(t)
# I should point out that the latter answer makes no sense since U(x) depends on s.
#Even using simplify on U(x) (as tomleslie does) doesn't change the behavior when done like this:
inttrans:-invlaplace(simplify(U(x)),s,t) assuming phi>0,x>0;
###################
Simpler example:
restart;
U(x):=x/(s^2+1);
inttrans:-invlaplace(U(x),s,t) assuming x>0 ;
                        U(x) Dirac(t)
%;
                        exp(sqrt(s))*Dirac(t)/(s^2+1)
inttrans:-invlaplace(U(x),s,t);
                        x*sin(t)
##The first part of following is OK
restart;
assume(x>0);
U(x):=x/(s^2+1):
inttrans:-invlaplace(U(x),s,t);
                            x~sin(t)
inttrans:-invlaplace(U(x),s,t) assuming x>0; #Again nonsense:
                           x~Dirac(t)
                           ----------
                              2      
                             s  + 1  

# When using assuming (as opposed to assume) I think the x in U(x) is temporaily given another name x~ (or some such thing) and U(x~) is not found in the remember table for U. Thus U(x~) is just an unknown quantity apparently not dependent on s, so the output will first be U(x~)*Dirac(t) but as it exits `assuming` x~ is replaced by x. 

#The comment above has been edited along the way.

@Markiyan Hirnyk I fail to see the relevance to the present situation in the link you give.

@Markiyan Hirnyk You may look up the Picard-Lindelõf Theorem in the excellent book by Philip Hartman, Ordinary Differential Equations, 1964, or any decent book past the most elementary level. Alternatively, take a look at
https://en.wikipedia.org/wiki/Picard%E2%80%93Lindel%C3%B6f_theorem

@Markiyan Hirnyk I shall make an attempt at an explanation.

The original problem is
ode1:=diff(y(t),t) = f(t,y(t),z(t)); # (1)
eq:=g(t,y(t),z(t)) = 0; # (2)
with y(0)=y0.
For mu>0 consider with h(t) = g(t,y(t),z(t)) the ode
ode3:=-mu*diff(h(t),t) = h(t); # (3)

Suppose that (1) and (3) with y(0)=y0 and z(0) = z0 has a solution on the interval [0,T), where T > 0. This is certainly the case under mild smothness assumptions. In some cases we may take T = infinity, but not necessarily (see example below).
We then have from ode3 (which is linear in h) that h(t) = h(0)*exp(-t/mu) for all t in the interval [0,T).
In the case where we may take T = infinity we see that h(t) -> 0 exponentially. For small values of mu h(t) is soon almost zero, so that the constraint equation eq is (almost) satisfied. Even if T < infinity we still have that h(t) decreases in absolute value. For a proper choice of mu it may become acceptably small fast enough.

Simple example of system where T cannot be taken to be infinity:
restart;
ode:=diff(y(t),t)=y(t)^2+z(t);
eq:=y(t)-z(t)=0;
dsolve({subs(z=y,ode),y(0)=1}); #Maximal T below is going to be near ln(2)
res:=dsolve({ode,eq,y(0)=1},numeric,method=rkf45_dae); #Just using dae-method in Maple
plots:-odeplot(res,[t,y(t)],0..0.69);
#Taking a somewhat small mu:
mu:=.01:
h:=lhs(eq);
ode2:=-mu*diff(h,t)=h;
res2:=dsolve({ode,ode2,y(0)=1,z(0)=0.8},numeric);
res2(.1); #Constraint (y=z) roughly satisfied at t = .1
plots:-odeplot(res2,[t,h],0..0.7);
h0:=eval(h,res2(0));
plots:-odeplot(res2,[t,h-h0*exp(-t/mu)],0..0.7);

@Markiyan Hirnyk I agree that there is no point (or difficult to see the point) in using a range like t=0..10^300.

Thanks for testing this in Mathematica. I wonder how NDSolve manages to escape the problem of x(t) becoming negative.
(Actually the same holds for Maple's dsolve with method=dver78 or taylorseries: see below):
It would be cheating if NDSolve actually ignored Abs, so that the equation became diff(x(t),t)=-x(t). For that equation there is no problem, since the equilibrium solution x(t) = 0 is asymptotically stable. The problem for the original equation diff(x(t),t)=-abs(x(t)) is that the equilibrium solution x(t) = 0 is unstable. More precisely stable from the right (x>0), but unstable from the left (x<0) in such a way that if x(t1) < 0 for some t1, then x(t) -> infinity for t -> infinity.

######
Briefly I tried the various methods in Maple (with default settings otherwise):

M:=[rkf45,ck45,rosenbrock,dverk78,lsode,gear,taylorseries]:
T:=100;
for m in M do
  try
    res:=dsolve({ode,x(0)=1},numeric,method=m);
    print(m,subs(res(T),x(t)));
  catch:
    print(lasterror)
  end try
end do;
#We see that for dverk78 and taylorseries x(t) remains positive when T=100.
The result for dverk78 is pretty good: Very close to evalf(exp(-100)).
Actually we can go higher for those two methods.

@maple fan Since the numerical integrator (with the default method = rkf45 and in Maple) starts getting negative values at time roughly t=34.4 the fact that dsolve doesn't want to continue beyond t=739.5991 is just because the value of x(t) is huge (and negative):
restart;
ode:=diff(x(t),t)=-abs(x(t));
res:=dsolve({ode,x(0)=1},numeric);
res(1000);
res(739.59941); #  x(t) = -1.35673646011395*10^304
## Using events for finding the "zero" of x(t):
resE:=dsolve({ode,x(0)=1,T(0)=0},numeric,discrete_variables=[T(t)::float],events=[[x(t)=0,T(t)=t]]);
resE(1000);
resE(770.28981);
###
Does Mathematica get into negative values just as Maple? In other words is the value found for t = 10^308 (or just before) huge and negative or is it incredibly small and positive? That I think is a more interesting question than when either system decides to call it quits.

In your pdf-file you mention in the beginning of Example 3 that Maple has problems with that equation.
However, as mentioned by acer in the link you give to MaplePrimes in your post above:
http://www.mapleprimes.com/questions/203096-Solving-Nonlinear-ODE#answer211682

dsolve handles the equation if the method is specified as (e.g.) rkf45_dae and a consistent initial condition is given for D(y1)(0):

eq1:=diff(y1(t),t)^2+diff(y1(t),t)*(y1(t)+1)+y1(t)=cos(diff(y1(t),t));
subs(diff(y1(t),t)=dydt0,y1(t)=0,eq1);
z0:=fsolve(%,dydt0);
solx:=dsolve({eq1,y1(0)=0,D(y1)(0)=z0},numeric,method=rkf45_dae);

I'm wondering if you are tacitly assuming that your system is autonomous, or that at least the algebraic equations are autonomous, i.e. don't have any explicit t-dependence.
I'm thinking about the initial phase where the switch is 0 so y is constant which supposedly gives z time to find a value consistent with the y-value. This appears to me to be problematic if g(t,y,z) actually depends explicitly on t.
Your examples are all autonomous.

@Markiyan Hirnyk The word 'bifurcate' has a meaning in ordinary nontechnical language. It means "to cause to divide into two branches or parts", see e.g. http://www.merriam-webster.com/dictionary/bifurcate
What happens for the simple system considered for t=Pi/2 where y=1 and z=0 is that locally the solution y(t)=sin(t), z(t)=cos(t) can be continued in two ways on any interval [Pi/2,Pi/2+epsilon], epsilon > 0. Either as the same analytic expressions or as y(t)=1, z(t)=0.
In the ordinary language use of the word 'bifurcate' it is indeed a bifurcation: A dividing into two branches.

Above I used the word 'locally'. As my previously given example illustrated there are actually infinitely many solutions passing through (y,z) = (1,0) (and also with the restriction that y(0)=0, z(0)=1).

@krw2015 I tried the lines

with(DifferentialGeometry); with(Tensor); with(Tools);

DGsetup([t, r, theta, phi], M1, verbose);

in Maple 17.02 and in Maple 2015.1. I didn't receive any errors in either version.

@ It could be pointed out that if we don't restrict ourselves to the interval -Pi/2..Pi/2 but consider solutions valid for all t then the problem
diff(y(t),t)=z(t), y(t)^2+z(t)^2=1 with y(0)=0,z(0)=1
has infinitely many solutions with y continuously differentiable and z continuous.
An example:

restart;
ode:=diff(y(t),t)=z(t);
eq:=y(t)^2+z(t)^2=1;
ics:={y(0)=0,z(0)=1};
dsolve({ode,eq});
Y:=piecewise(t<Pi/2,sin(t),t<4,1,t<4+Pi,sin(t+Pi/2-4),-1);
Z:=diff(Y,t);
plot(Y,t=0..10);
plot(Z,t=0..10,thickness=3);

#########Added note ########
It is interesting to see what dsolve does to the problem if we solve eq for z(t) and choose the nonnegative solution:
ode2:=diff(y(t),t)=sqrt(1-y(t)^2);
dsolve({ode2,y(0)=0});
#We get y(t) = sin(t), which is correct on the interval -Pi/2..Pi/2, but not on any larger interval containing t=0.
#Now numeric solution:
res:=dsolve({ode2,y(0)=0},numeric);
plots:-odeplot(res,[t,y(t)],-5..5); #The solution is constant outside the interval -Pi/2..Pi/2 as one would expect.

@Arthurok What matters is that f1(y) depends on y only and that f2(x) depends on x only.
The equations eq1 and eq2 gives you the only possible solutions for u(x,y) and v(x,y).

Inserting those into pde3 and differentiating the resulting equation w.r.t x and y as I did above you get the equation
-(2/3)*(-72*nu*y+336*y)/h^3+48*y*nu/h^3 = (48*(1+nu))*y/h^3
which must be satisfied for (at least) a range of y-values. Notice that f1(y) and f2(x) are gone.

Now the difference between the right and left hand sides of this equation should then be identically zero (on an interval at least). But the difference is (as found above)
-16*y*(-17+3*nu)/h^3

This is identically zero if and only if nu = 17/3. So in general the problem has no solution.

If on the other hand nu=17/3 then pdsolve can find the solutions:
res:=pdsolve(eval({pde1,pde2,pde3},nu=17/3));

However, your initial value problem has no solution:
pdsolve(eval({pde1,pde2,pde3,u(L,0) = 0,v(L,0) = 0,D[1](u)(L,0) = 0},nu=17/3));

Notice that the * you had in the third condition makes no sense. The derivative of u w.r.t. the first variable (x) at (x,y) = (L,0) is D[1](u)(L,0).

@Rouben Rostamian  I would like to call attention to the fact that I have edited the code for the two events. The changes are described in the edited version. One is logical the other is numerical.

I just tried copy and paste on your list of lists and replaced the comma at the very end with a semicolon.
I got one execution gropup, but as you say several prompts. Those latter prompts don't matter. It executed fine in Maple2015, worksheet mode and 1d input.

@Rouben Rostamian  First of all thank you very much for clearing this up.
I think I have a solution now using 'events'.
For convenience in changing parameters (e.g. g or P) I also use the parameters option, but that is for convenience only and can be left out. If so parameters must be assigned ahead of using dsolve.
#Edited relative to the original version:
1. Changed the second event from [a/mu/g=-1..1,[b(t)=1,n(t)=signum(v(t))]] to [a/mu/g=-1..1,[b(t)=1,n(t)=signum(a)]], because when abs(a)>mu*g the velocity will increase from its zero value if a>0 and decrease from zero if a < 0. Thus n(t) will have the value signum(v(t)) right after the event. Logically this makes sense. The success of the original could be called accidental.
2. Changed the first event from [v(t)=0,[b(t)=0,n(t)=0]] to [v(t)=0,[b(t)=0,n(t)=0,v(t)=0]] to ensure a perfect zero for v(t). This is done for numerical reasons.


restart;
X := t -> P*cos(omega*t);   # P for amplitude
#Let V := dX/dt, A := dV/dt, v := dx/dt, a := dv/dt.
#The equations of motion are de1 and de2:
de1 := diff(x(t),t) = v(t);
A:=diff(X(t),t,t);
a:=-mu*g*n(t)-A; #For n(t) see the dsolve command
de2:=diff(v(t),t)=b(t)*a; #For b(t) see the dsolve command
##
res:=dsolve(
       {de1,de2,x(0)=0.5,v(0)=0.25,b(0)=1,n(0)=1},numeric,
       discrete_variables=[b(t)::boolean,n(t)::integer],
       events=[[v(t)=0,[b(t)=0,n(t)=0,v(t)=0]],[a/mu/g=-1..1,[b(t)=1,n(t)=signum(a)]]],
       parameters=[P,omega,mu,g]
       );
res(parameters=[P=1,omega=1,mu=0.8,g=1]); #setting parameters
plots:-odeplot(res,[t,v(t)],0..15,thickness=3);
plots:-odeplot(res,[t,x(t)],0..15,thickness=3);
plots:-display(Array([plots:-odeplot(res,[t,n(t)],0..15,thickness=3),plots:-odeplot(res,[t,b(t)],0..15,thickness=3)]));
##The graph for v(t) with the parameters set as above:



###########################
#With parameters = [P=1,omega=1,mu=0.8,g=9.81] we get what you described: The house comes to relative rest and remains at relative rest.
#Now changing also P:
res(parameters=[P=10,omega=1,mu=0.8,g=9.81]);
plots:-odeplot(res,[t,v(t)],0..15,thickness=3);
plots:-odeplot(res,[t,x(t)],0..15,thickness=3);
plots:-display(Array([plots:-odeplot(res,[t,n(t)],0..15,thickness=3),plots:-odeplot(res,[t,b(t)],0..15,thickness=3)]));