As already mentioned by Markiyan this output is due to numerical roundoff.
The exact solution is x(t) = exp(-t), which obviously never becomes zero.

Your equation
ode:=diff(x(t),t)=-abs(x(t));
reduces to diff(x(t),t) = -x(t) on any interval on which x(t) > 0. Since x(0) = 1 > 0, the solution will be x(t) = exp(-t) as long as exp(-t) > 0, which it is for all t.

It is interesting that dsolve (in Maple 2015) returns 2 solutions, but one of those (x(t) = exp(t) ) is obviously wrong!

dsolve({ode,x(0)=1});
combine([%]);
                [x(t) = exp(t), x(t) = exp(-t)]

I checked in Maple 12. It gave only the correct solution.

Leaving out the initial condition in Maple 2015 gives
dsolve(ode);
                      exp(-t)                   
               x(t) = -------, x(t) = exp(t) _C1
                        _C1                     
This is only OK if interpreted correctly: the first _C1 is positive and the second _C1 is negative (or zero). Thus the C's should have had assumptions as are used in outputs from solve.

Doing the same in Maple 12 returns a much better result on implicit form and also better when using the option 'explicit'.
In Maple 2015 one could try using the option 'implicit':
dsolve(ode,implicit); #As in Maple 12
dsolve({ode,x(0)=1},implicit); #Using x(0)=1, but not entirely! The piecewise should have been simplified.

I shall submit an SCR about dsolve({ode,x(0)=1}).

Take a look at
?LinearAlgebra[RowOperation]

if you haven't done so already.

Try this simple example and take a look at the output (besides the warning):

p:=proc(x) y:=x+1 end proc;

#You notice that the output is actually the following, where y is declared local, something you ought to have done yourself. If you actually want y to be a global variable, then replace the word 'local, by 'global'.

p:=proc(x) local y; y:=x+1 end proc;

The system is inconsistent:

restart;
pde1:=diff(u(x, y), x) =-(2/3)*(3*h^3*nu+9*h^2*nu*y-12*nu*y^3+36*x^2*y+56*y^3)/h^3;
pde2:=diff(v(x, y), y) = (2/3)*(36*nu*x^2*y+56*nu*y^3+3*h^3+9*h^2*y-12*y^3)/h^3;
pde3:=diff(u(x, y), y)+diff(v(x, y), x) =-(6*(1+nu))*x*(h^2-4*y^2)/h^3;
pdsolve({pde1,pde2,pde3}); #Notice the warning
#Now by hand:
eq1:=map(int,pde1,x)+(0=f1(y));
eq2:=map(int,pde2,y)+(0=f2(x));
eval(pde3,{eq1,eq2});
diff(%,x,y);
simplify((rhs-lhs)(%)); #Not identically zero

Since you don't give the initial and boundary conditions the following will in part be an example.

restart;
pde:=diff(c(x,t),t)=diff(c(x,t),x,x)-diff(c(x,t),x)-c(x,t);
pdsolve(pde); #Variables separated
pdsolve(pde,build); #Variables separated and the odes solved
pdsolve({pde,c(0,t)=0,c(1,t)=0}); #Example boundary conditions
evalc(%);
about(_Z1); #_Z1 may be some other _Zn if you repeat the code without a restart
#Thus the solution will be of the form
C:=exp(x/2-5*t/4)*Sum(b(n)*sin(n*Pi*x)*exp(-n^2*Pi^2*t),n=1..infinity);
#Now the intial condition c(x,0)=f(x) (whatever it is) must be met:
f(x)=eval(C,t=0);
##
## The solutions give us the idea of making a simple change of dependent variable, which turns 'pde' into the heat equation:
PDEtools:-dchange({c(x,t)=u(x,t)*exp(x/2-5*t/4)},pde,[u]);
op(solve(%,{diff(u(x,t),t)}));
## Thus you can use tricks and methods from textbook examples.

You could use piecewise:

restart;
k := 10; m := 1; g := 10; mu := .2;
Xr := U*t-x(t);
ode1 := m*diff(x(t), t, t) = k*Xr-piecewise(mu*m*g<k*Xr,mu*m*g,0);
U:=1:
res:=dsolve({ode1,x(0)=0,D(x)(0)=0},numeric);
plots:-odeplot(res,[t,x(t)],0..15);
plots:-odeplot(res,[t,diff(x(t),t)],0..15);

The right hand side of ode2 can be written as -mu*m*g*signum(v)-m*A*omega^2*cos(omega*t).
That expression is discontinuous in v at v=0. That becomes a problem when v(t) hits zero (which it does with your initial conditions). When v(t) hits zero the expression changes with the amount 2*mu*m*g.

You could try using dsolve with events:

odeE:=diff(v(t),t)=-0.8*9.8*(2*b(t)-1)-cos(t)^2;
res:=dsolve({odeE,ode1,v(0)=0.25,x(0)=0.5,b(0)=1},numeric,discrete_variables=[b(t)::boolean],events=[[v(t)=0,b(t)=0]]);
plots:-odeplot(res,[t,v(t)],0..1,thickness=3);
plots:-odeplot(res,[t,x(t)],0..1,thickness=3);

The introduction of a floating point number like 1.0 is contagious. Try:
Pi*1.0;
#You will see that you don't get Pi and you don't get its infinite (!) decimal expansion either (obviously).
# sin(Pi*1.0) is going to give you the same as
Pi*1.0;
sin(%);
##To get 0 (well in fact the float 0.) you can do
evalf[11](sin(Pi*1.0));
# assuming that you are using Digits=10 (the default). Any number higher than 11 will do too.

Here is an example using only 3 points. You can just add 3 more points.
Before the animation command I give an image at time t=0:

restart;
with(plots):
f1:=cos(t)/4: f2:=sin(t)/5: f3:=sin(2*t+1)/6:
pts:=[[0,0,f1],[1+f2,1,1],[-1,0,f3]];
ptsL:=[pts[],pts[1]];
pointplot3d(eval(pts,t=0),symbol=solidsphere,symbolsize=50); p1:=%:
pointplot3d(eval(ptsL,t=0),connect,thickness=3,color=black); p2:=%:
display(p1,p2);
#Now animating from t = 0 to t = 50:
animate(display,['pointplot3d'(pts,symbol=solidsphere,symbolsize=50),'pointplot3d'(ptsL,connect,thickness=3,color=black)],t=0..50,frames=100);
##Notice the unevaluation quotes around pointplot3d in the animate command. They are there to prevent premature evaluation, i.e. that pointplot3d sees pts and ptsL as containing a symbol t, thus it cannot plot anything.

Assuming that you want real solutions (of which there are infinitely many in all 3 cases) you can use fsolve as follows.
Notice that because cos and cosh are both even, and tan and tanh are both odd, we can restrict ourselves to the positive real axis.

restart;
eq1:=cos(x)*cosh(x)=1;
plot((lhs-rhs)(eq1),x=0..15,-1..1);
fsolve(eq1,x=4..5); #Finds a root between 4 and 5
fsolve(eq1,x=14..15);
eq2:=cos(x)*cosh(x)=-1;
plot((lhs-rhs)(eq2),x=0..15,-1..1);
fsolve(eq2,x=10..12);
eq3:=tan(x)=tanh(x);
plot((lhs-rhs)(eq3),x=0..15,-1..1,discont=true);
fsolve(eq3,x=9..11);
########### solve ########
If you try this version of solve:
res:=solve(eq1,x,AllSolutions);
##
Then you will see that the result is a RootOf-expression:
RootOf(_Z-(2*I)*Pi*_Z5+2*arccosh(1/cos(_Z))*_B4-arccosh(1/cos(_Z)))
You can say that eq1 has just been rewritten some.
You will never get closer to an actual formula for all the roots.
If we use other notation we can write the argument to RootOf as
expr:=x-(2*I)*Pi*n+2*arccosh(1/cos(x))*b-arccosh(1/cos(x));
## _Z5 (or a different number than 5) stands for an arbitrary integer. _B4 can be 0 or 1.
## I have used n and b for those two. Instead of the unknown _Z I have used the name you use in eq1, x.
plot(eval(expr,{n=0,b=0}),x=0..15); #An example
fsolve(eval(expr,{n=0,b=0}),x=4..5); #Compare with the root found earlier

Looking at
indets(C21,name);
after running your worksheet P2.mw I found that theta occurs as just that and also as theta[1] and theta[2].
Furthermore r occurs as just r, but also as r[x], r[y], and r[z].

In general this is a bad idea in Maple because of the special meaning attached to (e.g.) theta[1].

I tried changing theta[1] to theta1, r[x] to rx, and the other 3 similarly.
Not being a friend of 2d-input I copied the whole thing to a new worksheet and converted to 1d input. There I could make the changes with Find/Replace in the Edit menu. (If anybody else is going to try that remember also to replace in the same way semicolon by colon, since the outputs generated are quite large and may force you to shut Maple down).
I didn't receive any errors when running the entire new worksheet.

I asked fsolve to solve for cos and sin as follows:

restart;
Digits:=30:
f1 := -8100+(-30+70*cos(t1)-40*cos(t2))^2+(-70*sin(t1)+40*sin(t2))^2;
f2 := (-20-80*cos(t3))^2+(-15+70*cos(t1)+10*cos(t1+t))^2+(-70*sin(t1)-10*sin(t1+t)+80*sin(t3))^2-5625;
f3 := (-20-80*cos(t3))^2+(15+40*cos(t2)+10*cos(t1+t))^2+(-40*sin(t2)-10*sin(t1+t)+80*sin(t3))^2-5625;
f4 := 10*cos(t1+t)*(30-70*cos(t1)+40*cos(t2))-10*sin(t1+t)*(70*sin(t1)-40*sin(t2));

idts:=indets({f1,f2,f3,f4},function);
S:=idts=~{c1,c2,c3,c4,s1,s2,s3,s4};
eqs:=subs(S,{f1,f2,f3,f4});
eqs2:={seq(cat(s,k)^2+cat(c,k)^2=1,k=1..4)};
r1:=fsolve(eqs union eqs2);
#Now that we know cos and sin for t1,t2,t3, and t4=t+t1 we can find the infinitely many solutions for (t1,t2,t3,t4).
#I didn't find anymore by using avoid:
fsolve(eqs union eqs2,{c1,c2,c3,c4,s1,s2,s3,s4},avoid=r1);
############Trying also solve which finds other solutions:
resS:=evalf(solve(eqs union eqs2));
subs(resS[1],eqs union eqs2); #Seems OK
subs(resS[2],eqs union eqs2); #Seems OK
subs(r1,eqs union eqs2); #Seems OK
############# Many more:
allvalues([solve(eqs union eqs2)]):
res:=evalf(%):
###The real ones:
resR:=remove(has,[res],I);
map(nops,resR); #resR is a list containing 8 lists each consisting of 2 sets of solutions
resR[1]; #Example
resR[1,1]; #Example
#Checking the solutions:
N:={}:
for i from 1 to 8 do
  for j from 1 to 2 do
    N:=N union subs(resR[i,j],eqs union (lhs-rhs)~(eqs2))
  end do
end do;
N;
fnormal(N,21);

That you get the empty set is not surprising since (x-2)/(x-2) automatically simplifies to 1 before any actual evaluation.
This is evidenced by trying unevaluation quotes around the input:
'(x-2)/(x-2)'; #You get 1.
#The package InertForm was introduced in Maple 18:
with(InertForm);
u:=InertForm:-Parse("(x-2)/(x-2)");
InertForm:-Display(u);
indets(u,name);
Value(u);

This is clearly difficult to deal with.
I changed blt from 15.5 to just 3. I had only success with s = 0.
It appears that theta for all practical purposes is zero for eta > 0.003. Thus I tried solving Eq1 with theta = 0.
That was very easy. After that my intention was to use the results for f to solve Eq2.

restart;
Digits:=15:
with(plots):
#blt:=15.5:
A:=1:
M:=1:
lambda:=1:
fw:=1:
rhos:=8933:
rhof:=997.1:
rhob:=14918.11:
rhob2:=20939.1:
kn:=0.613:
kf:=401:
cps:=385:
cpf:=4179:
Pr:=6.2:
Eq1 := (1/(1-s)^(2.5))*diff(f(eta),eta,eta,eta)-(1-s+s*(rhos/rhof))*(diff(f(eta),eta)^(2)-f(eta)*diff(f(eta),eta,eta)+A*(diff(f(eta),eta)+(1/2)*eta*diff(f(eta),eta,eta)))-M*diff(f(eta),eta)+(1-s+s*((rhob/rhob2)))*lambda*theta(eta)=0;
Eq2 :=(1/Pr)*((kn/kf)/(1-s+s*(rhos*cps/rhof*cpf)))*diff(theta(eta),eta,eta)+f(eta)*diff(theta(eta),eta)-diff(f(eta),eta)*theta(eta)-A*(2*theta(eta)+0.5*eta*diff(theta(eta),eta))=0;
blt:=3:
bcs1 := f(0) = fw, (D(f))(0) = 1, (D(f))(blt) = 0, theta(0) = 1, theta(blt) = 0;
L := [0,0.05,0.1,0.2];
R:=dsolve(eval({Eq1, Eq2, bcs1}, s = 0),numeric, output = listprocedure,initmesh=1024,maxmesh=8192);
odeplot(R,[eta,theta(eta)],0..0.003);
odeplot(R,[[eta,f(eta)],[eta,theta(eta)]],0..blt);
#################
##Now the approximation obtained by setting theta = 0 in Eq1.
## I'm saving the procedures for f, f', f'' for use later:
blt:=8: #I could get away with blt=8
AP:=NULL:
for k from 1 to 4 do
  Rf[k] := dsolve(eval({Eq1f, f(0)=1,D(f)(0)=1,D(f)(blt)=0}, s = L[k]),numeric, output = listprocedure,AP);
  AP:=approxsoln=Rf[k];
  F[k],F1[k],F2[k]:=op(subs(Rf[k],[seq(diff(f(eta),[eta$i]),i=0..2)]))
  end do;
display(seq(odeplot(Rf[k],[eta,f(eta)],0..blt),k=1..4),insequence);
################
#Now I tried the following loop, but had success only with k=1, i.e. s=0:
AP:=NULL:
for k from 1 to 4 do
 Rt[k]:=dsolve({subs(diff(f(eta),eta)=F1[k](eta),f(eta)=F[k](eta),eval(Eq2,s=L[k])),theta(0)=1,theta(0.003)=0},numeric,known=[F[k],F1[k]],AP);
AP:=approxsoln=Rt[k]
end do;
odeplot(Rt[1],[eta,theta(eta)],0..0.003);
###############FINALLY I tried an approximation on Eq2 which uses that theta decreases steeply to zero.
## I'm using that Eq2 can be approximated by K*theta''+f*theta' = 0, but also by
## K*theta''+(f*theta)' = 0 since f '*theta is relatively small.
## Integrating once using that theta and theta' are zero for eta = 0.003 we get the first order ode:
Eq20 :=(1/Pr)*((kn/kf)/(1-s+s*(rhos*cps/rhof*cpf)))*diff(theta(eta),eta)+f(eta)*theta(eta)=0;
#This can be solve symbolically in terms of f:
solTH:=dsolve({Eq20,theta(0)=1},theta(eta));
#Now we can plot theta
for k from 1 to 4 do
  if k=1 then r:=0.003 else r:=1e-8 end if;
  pp[k]:=plot(eval(rhs(solTH),{s=L[k],f=F[k]}),eta=0..r,thickness=5)
end do:
display(Array([seq(pp[i],i=1..4)]));
###In fact we may as well replace f with 1 in solTH. It makes no difference of any importance to the results, but has the grat advantage that solTH then can be integrated without knowing f:
solTH2:=value(eval(solTH,f=1));
#Now on the other hand we could replace theta in Eq1 with that explicit expression:
Eq1f2:=subs(solTH2,Eq1);
##Now use a dsolve loop as above taking blt=5:
blt:=5:
AP:=NULL:
for k from 1 to 4 do
  Rf2[k] := dsolve(eval({Eq1f2, f(0)=1,D(f)(0)=1,D(f)(blt)=0}, s = L[k]),numeric, output = listprocedure,AP,maxmesh=8192,abserr=1e-5);
  AP:=approxsoln=Rf2[k];
  #F[k],F1[k],F2[k]:=op(subs(Rf2[k],[seq(diff(f(eta),[eta$i]),i=0..2)]))
  end do;
The graphs are really indistinguishable from the ones obtained using Eq1f with blt=8 (ignoring eta>5).

My conclusion would be to use solTH2 as the solutions for theta and as the solutions for f the ones obtained from Eq1f.

I don't think you are going to find a solution with the values of R and delta that you have. So either R has to be smaller or delta has to be larger.
Instead of going through a loop solving the ode system 3000 times I used the parameters option in dsolve and also fsolve.
For convenience I bring the entire code.

restart;
R:=0.3:
theta_min:=Pi/6:
theta_max:=Pi/2:
betta_max:=evalf(Pi/180*80);
p:=2*10^5:
theta0:=s->Pi/3/s_end*s+Pi/6:
r0:=s->R*sin(theta0(s)):
s_end:=evalf(R*(theta_max-theta_min)):
sol1:=solve({sin(betta_max)=c/r0(0)},{c});
const1:=0.1477211630;
betta0:=s->arcsin(const1/r0(s)):
betta:=s->arcsin(r(s)/r0(s)*sin(betta0(s))):
A:=s->cos(betta(s))/cos(betta0(s)):
T1:=s->rT1(s)/r(s):
T2:=s->T1(s)*tan(betta(s))^2:
delta:=0.001:
sys := diff(rT1(s), s)-A(s)*T2(s)*cos(theta(s)),diff(theta(s), s)-A(s)/T1(s)*(p-T2(s)*sin(theta(s)/r(s))),diff(r(s),s)-A(s)*cos(theta(s)),diff(z(s),s)-A(s)*sin(theta(s));
rT1_n:=p*Pi*r_min^2/2/Pi/sin(theta_min);
##### No changes made above to the lines included. ############
## Now use the parameters option with r_min as parameter:
F := dsolve({sys,rT1(0)=rT1_n, theta(0)=theta_min,r(0) = r_min, z(0) = 0}, numeric,output=listprocedure,parameters=[r_min]);
r_ravn:=subs(F,r(s)); #The procedure for evaluating r(s) extracted
#Testing. Setting the parameter:
F(parameters=[0.152]);
F(s_end); #Values at s_end
##Plotting (scaling so all 3 appear) (still r_min=0.152)
plots:-odeplot(F,[[s,10*r(s)],[s,rT1(s)/20000],[s,theta(s)]],0..s_end);
## Another plot with separate plots (still r_min=0.152)
use plots in display(Array([seq(odeplot(F,[s,fu(s)],0..s_end),fu=[r,rT1,theta])])) end use;
## Making a procedure which sets the parameter and plots. Default all 3 functions as an array.
q:=proc(r_min,L::list:=[r,rT1,theta]) uses plots; F(parameters=[r_min]);
  display(Array([seq(odeplot(F,[s,fu(s)],0..s_end),fu=L)])) end proc;
q(.001); #Plotting all with r_min = 0.001
q(0.1,[r]); #Plotting only r(s) with r_min=0.1
plots:-animate(q,[rmin],rmin=0.01..0.152,frames=100); #Animation (all 3)
plots:-animate(q,[rmin,[r]],rmin=0.001..0.152,frames=100); #Animation (only r)
##It appears that r stays below R-delta for all s = 0..s_end
## Thus I replace abs(r_ravn(s_end)-R)=delta by r_ravn(s_end)-R+delta=0 below.
## In the procedure pp the default is to use the existing delta, but another can be given as a second argument:
pp:=proc(r_min,dd::numeric:=delta) F(parameters=[r_min]); r_ravn(s_end)-R+dd end proc;
#Testing:
pp(.152);
pp(.152,.01);
pp(.152,.02);
#Now use fsolve. I expected no success with the default delta and didn't have any:
fsolve(pp,0.01..0.152); returns unevaluated
fsolve(rm->pp(rm,.02),0.01..0.152); #Success
### Finally to see that R-delta is too large we can plot r_ravn(s_end) as a funcion of r_min.
## The plot clearly shows that the maximal value for r_ravn(s_end) is for r_min=0.152:
vv:=proc(r_min) F(parameters=[r_min]); r_ravn(s_end) end proc;
plot(vv,0.01..0.152);
vv(0.152);