I think you are demanding too much in asking that U be 0 for t<=0 and 1 for t>=1. As you will see the requirement that U(0)=U0 already determines U (I call the new function P instead of U).

With eq3 and sys as you define them I did this:

L:=convert(indets(eq3,piecewise),list);
S:=L=~[diff(P(t),t),P(t)]; #Make sure of the order of the elements in indets
numer(subs(S,eq3));
EQ3:=isolate(%=0,diff(P(t),t)); #New eq3, but in P(t)
SYS:=subs(S,{sys[2..3]}); #New equations for x1 and x2
NewSys:=SYS union {EQ3}; #The new system
res:=dsolve(NewSys union {x1(0)=x10,x2(0)=x20,P(0)=U0},numeric);
plots:-odeplot(res,[[t,P(t)],[t,x1(t)],[t,x2(t)]],-3..8,view=0..1);

#####
##The new system is
{diff(P(t), t) = (-P(t)^2*x1(t)+2*P(t)*x1(t)*sqrt(x2(t))+P(t)^2*x2(t)-2*P(t)*x2(t)*sqrt(x1(t)))/(-2*x1(t)^(3/2)*P(t)+4*x1(t)^(3/2)*sqrt(x2(t))+2*x2(t)^(3/2)*P(t)-4*x2(t)^(3/2)*sqrt(x1(t))), diff(x1(t), t) = P(t)-sqrt(x1(t)), diff(x2(t), t) = P(t)-sqrt(x2(t))};

Since the allowable input type to the discont option is either boolean or list there is a whole lot of nonsensical input (agreeing with Carl!), which will pass the initial typecheck.
Somewhat surprising is it that it doesn't run into other difficulties.

type(a<=7,boolean);
plot(1/x, x=-1..1, discont= (a<=7) ,thickness=2);
plot(1/x, x=-1..1, discont=["Anything",7,9,13] ,thickness=2);
###
#It may be that somewhere the following logic is applied:

p:=proc(x::{boolean,list}) if x=false then fff else ttt end if end proc;
#Actually I think this is a good and economical way to handle the option discont. The harm done is not great: The plot is done.

#All but the first of these return ttt:
p(false);
p(FAIL);
p(true);
p(a<=7);
p(["Anything",7,9,13]);
p([]);

Here is a simple version using a do loop.
p:=proc(n) local k,s;
   s:=0;
   for k from 0 to n do
      s:=evalf(s+4*(-1)^k/(2*k+1))
   end do;
   s
end proc;

I don't understand what the term Gradient(1-f * gradient(g)) is supposed to mean.
But for purposes of showing some syntax it probably doesn't matter so much which interpretation I now choose.
In the following I shall assume (most likely incorrectly) that it was intended as
Gradient(1-f ). Gradient(g), i.e. the inner product of the two gradients.

To avoid introducing conjugates I write Gradient(1-f )^%T. Gradient(g) instead. Here A^%T means the transposed of A.
In defining the equations I make use of Gradient and Laplacian from the VectorCalculus package and also the construction 'use ... end use' for convenience:

use grad=VectorCalculus:-Gradient,L=VectorCalculus:-Laplacian in
  eq1:=diff(f(t,x,y,z),t)=f(t,x,y,z)*(1-f(t,x,y,z))-f(t,x,y,z)*h(t,x,y,z);
  eq2:=diff(g(t,x,y,z),t)=g(t,x,y,z)*(1-g(t,x,y,z))+grad(1-f(t,x,y,z),[x,y,z])^%T.grad(g(t,x,y,z),[x,y,z]);
  eq3:=diff(h(t,x,y,z),t)=g(t,x,y,z)-h(t,x,y,z) +L(h(t,x,y,z),[x,y,z]);
end use;

The following is syntactically correct but WARNING: do not attempt it. I had to stop the proces and restart my computer.

pdsolve({eq1,eq2,eq3});

That would be an attempt to solve analytically. No wonder it failed.

As for solving numerically, you are out of luck. Maple can only solve time-dependent pdes in two variables, see
?pdsolve,numeric

Z is a set. After the assignments to h3, x3, g2 two of the equations in Z are the same, thus Z now has only 4 elements. Thus numelems(Z) = 4.
The interesting thing is that when you reference Z[1] and Z[4] they are both 0 = 0.
This must mean that the Z in Z[i] doesn't get evaluated fully first and then the ith element taken, but that the full evaluation takes place after the ith element has been selected.
If you replace numelems in the second loop with 5, you will see all.
If you use a list instead there is no problem.

Simple example.
restart;
A:={a,b,c};
b:=a;
A1:=A;
seq(A[i],i=1..3);
seq(A1[i],i=1..3);
seq(A1[i],i=1..2);

ArrayInterpolation wants floating point values it says. So you just make sure it gets such.
Using evalf on the data will do it.
Don't use with( package ) inside a procedure.

InterProc := proc ()
local i; global A,variable1, variable2;
A := Matrix(1 .. 4);
variable1 := evalf([0, 100]);
variable2 := evalf([12, 20]);
for i from 1 to 10 do
 A(i, 1) := 3+i;
 A(i, 2) := 2*i+A(i, 1);
 A(i, 3) := A(i, 2)-1;
 A(i, 4) := CurveFitting:-ArrayInterpolation(variable1, variable2, A(i, 1));
end do;
end proc;

InterProc();


#To me it seems better to have the data be given as arguments and have A local:

InterProc := proc(variable1::list(float), variable2::list(float))
local i, A;
A := Matrix(1 .. 4);
for i from 1 to 10 do
 A(i, 1) := 3+i;
 A(i, 2) := 2*i+A(i, 1);
 A(i, 3) := A(i, 2)-1;
 A(i, 4) := CurveFitting:-ArrayInterpolation(variable1, variable2, A(i, 1));
end do
end proc;
A:=InterProc([0., 100.],[12., 20.]); #global A here (or another name)

I'm somewhat surprised that you didn't link to your previous question where I answered this very same question.
http://www.mapleprimes.com/questions/200857-Compact-Solution-To-An-Ode
To repeat:

I think you miscalculated some.
But the following does come up with a relatively compact answer. The boundary conditions for U are assuming that L > 0.

restart;
ODE:=(diff(T(x), x, x))+P*(S+a*(1-exp(-L*x))/L)*(diff(T(x), x))=0;
bcs:=T(0)=1,T(infinit)=0;
var:=z=P/L^2*exp(-L*x);
var2:=x=solve(var,x);
PDEtools:-dchange({var2,T(x)=U(z)},ODE,[z,U(z)]);
collect(%,diff,factor);
ode1:=isolate(%,diff(U(z),z,z));
dsolve({ode1,U(P/L^2)=1,U(0)=0});

simplex is an old package. In display it uses some old objects like matrix (not Matrix) and even a function MATRIX. Clearly done for display purposes mainly, but you can extract the parts.

with(simplex);
res:=display( { x + 3*y + z <= 0, w - 2*y - z  <= 2 } ); #Example from help page
A:=op(1,lhs(res)); #Extracting A
var:=op(2,lhs(res)); #Extracting the variable "vector"
b:=rhs(res); #The right hand side
type(A,matrix); #Yes!
type(A,Matrix); #No
whattype(var); #Surprise: function
op(0,var); #What function: MATRIX
op(1,var); #Argument to MATRIX
evalm(A)*var; #using the old evalm
type(evalm(b),matrix);#Yes
type(simplex,table); #Yes, so an old package

Your answer using the first set of equations is correct though, and the second set returns the same.

res:=solve([sys1, sys2, sys3, sys4, sys5, sys6, sys7, sys8, sys9], {t1,t2,t3,t4,t5,t6,t7,t8,t9}); #{} used instead of []

eval(osys,res);
simplify(%);
simplify(%) assuming t4>0;
#However, for t4 < 0 it is wrong:
simplify(%%) assuming t4<0;

#The given result is a special case:
eval(res,{t1 = 1,t2 = 0,t3 = 0,t4 = 1,t5 = 0,t6 = 1,t7 = 0,t8 = 0,t9 = 1});
#Results in the set {0=0,1=1}, i.e. identities.

One of the problems comes up right at the top.
You want x to be of type numeric, OK. But then when you call it as in proc_r(7) the first thing that happens is that the formal parameter x is replaced by 7.
Thus you want to make the assignment 7:=x0, which is clearly not allowed.

proc_r:= proc(x :: numeric)
local x0;
x := x0;
end proc;

proc_r(7);
You probably meant x0:=x instead.

After that is corrected you have another problem: break can only be used in a do loop.
You have break(s) outside as in plain and simple:

if 1 <= 2 then  break end if;

which throws an error.

S[1] is initialized as 0 by S:=Array(0..(x-1));

In this part of the code S[1] is not updated and the do loop runs (almost) twice since S[0] is not equal to T.
The second check of equality is between S[1] and T, i.e. between 0 and [3 0 0 ]
for y from 0 to j do  #do2
  print("ddd",S[y],T);
  if LinearAlgebra:-Equal(S[y],T) then
    found:=true;
    if (found=false) then
      S[j]:=T;
      print("ccc",S[j]);
      j:=j+1;
    end if:
  end if:
  end do; ###end do2

You may have meant an 'else' instead of the second if-statement?

The syntax is wrong: no interval  like x = 0..1 should be given.
Also as a first order ode you can only give one condition.

fx := diff(n(x), x)-a;
dsolve({fx,n(0) = 1});

First thing: spline is deprecated. You should use Spline as you probably also intended since you loaded the CurveFitting package. Then the last argument should be degree = 2 (if that is what you want).

Second: You didn't give an x - range as in plot(Piece2, x= 0..15). When you don't do that plot assumes that the range is x= -10..10. The view option you use only limits or expands what you see. It doesn't affect what points are computed.
Thus replace

plot(Piece2, view = [0 .. 15, -1 .. 1]);

with
plot(Piece2, x=0..15);

and
plot(Piece1, view = [0 .. 15, -1 .. 1]);
with
plot(Piece1, x=0..10]);

Because of the factor cosh(t)*sin(2*t) the numbers 5+7*I, 5-7*I, -5+7*I, and -5-7*I are roots of the characteristic equation. Because of the factor t^7 they are roots of order 8.
So we have:
restart;
y:=3*t^7*cosh(5*t)*sin(2*t);
chareq:=((lambda-(5+2*I))*(lambda-(5-2*I))*(lambda-(-5+2*I))*(lambda-(-5-2*I)))^8;
ode:=DEtools[diffop2de](chareq,x(t),[lambda,t])=0;
odetest(x(t)=y,ode);

Why change time in order to go forwards? dsolve has nothing against going backwards:

N:=1.2345: T:=1:
res:=dsolve({diff(S(t), t) = -K(t)*S(t)/N, diff(K(t), t) = K(t)*S(t)/N,S(T)=10,K(T)=20},numeric);
plots:-odeplot(res,[[t,S(t)],[t,K(t)]],0..T);