@Axel Vogt Why is the given integral equation not
y(t) = 1-h*Int(JacobiTheta3(0,exp(-Pi^2*s))*y(t-s)^4,s = 0 .. t) ?
That seems to me to be implied by his recursion.
The discussion in the previous question posed by Al86 left me somewhat confused, so maybe I'm missing something.

@tomleslie I apologize for editing my comment many times within about one hour yesterday.
If you read the comment as it ended up yesterday after my last edit, I think the reason for the weird behavior is the following:
When using U(x):=uuu; #Where uuu is some expression, e.g. x/(1+s^2).
the entry U(x)=uuu is written in the memory table of U.
When using
inttrans:-invlaplace(U(x),s,t) assuming x>0;
the first that happens (before evaluation of U(x)) is that x is replaced by another (local) name (think of it as x~).
After that evaluation of U(x~) takes place. First the memory table is searched, the index x~ is not found. Thus since no (general) definition of U has been done, U(x~) is just itself, i.e. evaluates to U(x~).
After that invlaplace gets to work resulting in U(x~)*Dirac(t).
Just before exiting `assuming` the local name x~ is replaced by the global x resulting in the output U(x)*Dirac(t).
Notice that if you evaluate that by doing just
%;
you get uuu*Diract(t).
###
So is this a bug? In my view it is not. It is due to the deliberate design of `assuming`.
As I also pointed out (and you mention) the mess could easily have been avoided by using U or Ux (or whatever name) instead of U(x).

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.