You can use one more eval:

restart;
f:=t->piecewise(t<0,0,t>=0 and t<z,t,t>z,z):
r:=convert(f(t),Heaviside):
r:=inttrans[laplace](r,t,s);
eval(r,z=0.5);
eval(%,laplace=inttrans[laplace]);
#or these:
subs({z=0.5,laplace=inttrans[laplace]},r);
%;
subs(z=0.5,laplace=inttrans[laplace],r);
%;
#But not this one, which causes an infinite recursion:
eval(r,{z=0.5,laplace=inttrans[laplace]});

I wonder why inttrans[laplace](g(t),t,s) when it returns unevaluated doesn't return inttrans[laplace](g(t),t,s) instead of laplace(g(t),t,s).

Here are two ideas you may consider:

sys:=diff(x(t),t)=sin(t)+sin(y(t)),diff(y(t),t)=cos(t*x(t)); #Example
res:=dsolve({sys,x(0)=0,y(0)=0},numeric,output=listprocedure);
plots:-odeplot(res,[[t,x(t)],[t,y(t)]],0..5);
X,Y:=op(subs(res,[x(t),y(t)])); #Numerical procedures to compute x(t) and y(t)
Optimization:-Maximize(X(t),t=0..5);
Optimization:-Maximize(Y(t),t=0..2);
#Be aware that Maximize only gives a local maximum, so you can also try:
Optimization:-NLPSolve(X(t),t=0..5,method=branchandbound,maximize);
Optimization:-NLPSolve(Y(t),t=0..5,method=branchandbound,maximize);

#Alternatively, just look for zeros of x' and y':
fsolve(eval(rhs(sys[1]),{x=X,y=Y}),t=4);
fsolve(eval(rhs(sys[2]),{x=X,y=Y}),t=1);

u:=0.569840290998053;
convert(u,rational);
convert(u,rational,exact);
?convert,rational

I suspect some confusion here.
Are the values of y only allowed to be real numbers in the interval 0..10?
If so then the fsolve command which has y in 0..10 is relevant although it only looks for roots of Re(t). But then the RootFinding:-Analytic is not because it searches for roots (of t) in the y-rectangle re=0..10, im=-1..10.
If y is required to be in the interval 0..10 then the roots of t are the same as the roots of BesselJ(0,y) since 
evalc(abs(t)) assuming y>0;
The first factor cannot be zero for y>=0. The zeros of t are then the zeros of BesselJ(0,y):
BesselJZeros(0.,1..5); #The first 5 zeros, only 3 of which are in 0..10.

To solve the boundary value problem you will have to solve numerically.
In order to do that you have to make sure that the highest derivatives of y1 and y2 can be isolated by solve. That is not the case. The solution to that is to differentiate the second equation w.r.t. x.

restart;
#Notice that y1 has to be y1(x):
dsys:={diff(y1(x),x$2)-A*x^2*y1(x)=B*diff(y2(x),x$3),diff(y2(x),x$2)-C*diff(y1(x),x)=F*x^2+G};
bcs:={y1(0)=0, (D@@2)(y2)(L/2)=0, D(y2)(0)=0, y2(L/2)=0};
dsolve(dsys); #You get an answer, but that involves integrals
solve(dsys,{diff(y2(x),x$3),diff(y1(x),x$2)}); Not possible
dsys[2]; #The second equation
eq1:=diff(dsys[2],x); #The new second equation
convert(dsys[2],D); #We now need an extra boundary condition for y1, but that can be found:
eval(%,x=L/2);
bcs1:=eval(%,bcs);
dsys2:=solve({dsys[1],eq1},{diff(y2(x),x$3),diff(y1(x),x$2)});
res:=dsolve(eval(dsys2 union bcs union {bcs1}, {A=1,B=2,C=1,F=3,G=4,L=2}),numeric);
plots:-odeplot(res,[[x,y1(x)],[x,y2(x)]],0..1);
##Comment: If you try solving the new system symbolically:
dsolve(dsys2);
you notice that you still have integrals.

Besides the integral problem mentioned in my comment you have another problem: The distribution Dirac (in fact its derivative) is in the equation PDE2.
Consider the following very simple examples.
ode:=diff(u(t),t)=Dirac(t-1);
dsolve({ode,u(0)=0});
res:=dsolve({ode,u(0)=0},numeric);
plots:-odeplot(res,[t,u(t)],0..2);
ode2:=diff(u(t),t,t)=Dirac(1,t-1);
dsolve({ode2,u(0)=0,D(u)(0)=0});
res2:=dsolve({ode2,u(0)=0,D(u)(0)=0},numeric);
plots:-odeplot(res2,[t,u(t)],0..2);

You will notice that the numerical algorithm will just treat Dirac as zero.
So I think you ought to reconsider the formulation of your problem.

The error message tells the story. It is also clearly stated in the help:
?pdsolve/numeric

Your problem obviously has the solution u(x,y,t) = exp(2*t+x+y), but surely you just wanted to test a numerical solving procedure?

Actually, pdsolve can solve your problem (but not numerically):
pdsolve({PDE1} union IC union BC);

I think you got numerical problems.
If you try
pde2:=isolate(pde,diff(u(x,y),x));
then the right hand side has two terms T(x)/u(x,y) and diff(u(x,y),y,y)*10^(-6)/u(x,y).
The first evaluates to -2.17*10^11 at x=0, the second to 5.6*10^(-8).
It is not surprising therefore that the algorithm has a problem getting started.
If you decrease the timestep (x is time) you get the algorithm working at least on a small interval:
sol := pdsolve(pde, {u(0, y) = 18, u(x, 0) = 0.6e-3, u(x, 0.25e-4) = 18}, numeric, timestep = 1e-12, spacestep = 1e-7);
sol:-plot3d(x = 0 .. 4e-11);

You could also try scaling time (x) and space (y).

Not unless you write the code yourself, see the question raised by Chia in
http://www.mapleprimes.com/questions/202107-Solution-Structure-Of-Pdsolve-With-BC

during a discussion on a bug in pdsolve/numeric (not relevant here), where I pointed out that in the help page for pdsolve/numeric it say that the method used is time based and that that must mean that
the algorithm "marches" forward in time (t) from known initial conditions all along satisfying the boundary conditions at each end of the x-interval. You get no complaints therefore if you impose initial and boundary conditions which allow that to happen.

The example used there was
pde:=diff(u(x,y),x,x)+diff(u(x,y),y,y)=0;
bcs:={u(x,0)=0,u(x,1)=0,u(0,y)=y,u(1,y)=y^2};
pdsolve(pde,bcs,numeric); #Error about not handling elliptic pdes
#and then with icbs allowing the time based solved to "work":
ibcs:={u(x,0)=0,D[2](u)(x,0)=0,u(0,y)=y,u(1,y)=y^2};
res:=pdsolve(pde,ibcs,numeric,timestep=.001,spacestep=.001);
Quoting also the comment:
Whether it is a reasonable problem (ill-posed, or interesting) is another matter.

You can use unapply:

f:=unapply(MY_EXPRESSION, seq(t[i],i=1..N) );

Simple example:
N:=4:
MY_EXPRESSION:=add(t[i]*i,i=1..N);
f:=unapply(MY_EXPRESSION, seq(t[i],i=1..N) );

F2 will give sum2 fractional values on occasion. But
7/2 mod 24;

returns an error.
An unrelated error: You probably meant deltalist[i]<>0 in F2 instead of just deltalist.

Without trying to understand what you are up to I tried F1,F2, and F3 on
L:=[seq(1..35)];
like this
F1(L,L);

Also I tried
myProc40(L);
which returned the error:

Error, (in F2) the modular inverse does not exist

You probably should tell us what your input to myProc40 was.

I tried to follow Kamke.

restart;
ode:=diff(y(x), x,x)=(a*y(x)^2+b*x*y(x)+c*x^2+alpha*y(x)+beta*x+gamma)^(-3/2);
#Kamke 6.13:
PDEtools:-dchange({y(x)=u(x)-(b*x+alpha)/2/a},ode,[u]);
expand(%);
#Kamke 6.101
ode2:=(A*x^2+B*x+C)^(3/2)*diff(u(x),x,x)=f(u(x)/sqrt(A*x^2+B*x+C));
f:=q->1/(a*q^2+1)^(3/2);
PDEtools:-dchange({u(x)=sqrt(A*x^2+B*x+C)*q(x)},ode2,[q(x)]);
2*diff(q(x),x)*%;
map(int,%,x);
#Result in Kamke 6.101
ode3:=collect(%,[diff,q]); #q' is squared, we really get two odes:
S:=solve(ode3,{diff(q(x),x)});
ode4:=op(1,S[1]);
ode5:=op(1,S[2]);
res4:=dsolve(ode4);
res5:=dsolve(ode5);
odetest(res4,ode4);
odetest(res5,ode5);
sol:=dsolve(ode,implicit) assuming a::NonZero; #The assumption (at least) seems to speed up the computation
nops({sol});
odetest(sol[1],ode); #Goes on for a long time. My patience ran out.
#Obviously by "unsolvable" is not meant that the ode has no solution, since the usual existence results to initial value problems apply to that one as well. Example:
res:=dsolve({eval(ode,{a=1,b=1,c=1,alpha=1,beta=1,gamma=1}),y(0)=0,D(y)(0)=1},numeric);
plots:-odeplot(res,[x,y(x)],-3..2); 

By adding the first two equations we see that n(x,t)*v(x,t) depends on t only:
restart;
sis := [diff(c(x, t)*v(x, t)*n(x, t), x)+diff(c(x, t)*n(x, t), t) = diff(n(x, t), t), diff((1-c(x, t))*n(x, t)*v(x, t), x)+diff((1-c(x, t))*n(x, t), t) = 0, diff(n(x, t), t) = (48*(1/1000))*n(x, t)*c(x, t)*(1-n(x, t))*v(x, t)];
sis[1]+sis[2];
simplify(%);
(lhs-rhs)(%)=0;
#We let f(t) = n(x,t)*v(x,t). Thus we really only have these equations:
eq1:=eval(sis[1],v(x,t)=f(t)/n(x,t));
eq2:=eval(sis[3],v(x,t)=f(t)/n(x,t));
#Isolating the time derivatives:
solve({eq1,eq2},{diff(n(x,t),t),diff(c(x,t),t)});
#Notice now that if c(x,0)=0, then so is diff(c(x,t),x) for t=0. Thus from the equation just found for diff(c(x,t),t) it follows that c(x,t) = 0 for all t>0 too. From the equation found for diff(n(x,t),t) it then follows that n(x,t) is constant, i.e. equal to 0.4 for all x and t.

You can use the textbook method of separation of variables after a simple change of variables.

restart;
pde:=diff(T(x, y), x)-k*diff(T(x, y), y, y) = 0;
bcs:=D[2](T)(x,0)=a,D[2](T)(x,c)=b;
ch:=T(x,y)=u(x,y)+a*y+(b-a)/c*y^2/2+k*(b-a)/c*x;
PDEtools:-dchange(ch,pde,[u]);
simplify(%);
#Thus u satisfies the same pde as T does, but has the advantage of satisfying homogeneous boundary conditions:
diff(ch,y);
convert(%,D);
eval(%,y=0); #Thus D[2](u)(x,0)=0
eval(%%,y=c);# Thus D[2](u)(x,c)=0.
##Now follow any textbook in separating variables:
pdsolve(pde2);
#Solve the ode for _F2(y). Impose homogeneous boundary conditions on the derivatives of _F2. Those can only be satisfied for _c1<=0. Thus your solutions for _F2 are cosines for _c1<0 and a polynomial of degree 1 for _c1=0. Any linear combination of such terms still satisfies the boundary conditions.
Now you need initial conditions to proceed.

If you do have a concrete initial condition of the form T(0,y)=g(y) you can also solve numerically.

Is this what you have in mind?
I use z as the one name for  zeta and c. Not knowing gDiffusion or KDiffusion I set them at 1.

pde:=diff(z(x, t), t) = piecewise(x<2000,KDiffusion*10^6,gDiffusion*10^5)*diff(z(x, t), x,x);
ibc := {z(0, t) = .4, z(x, 0) =piecewise(x<2000,.4,0), D[1](z)(3000, t) = sin(1/100*t)};
sol:=pdsolve(eval(pde,{gDiffusion=1,KDiffusion=1}),ibc,numeric,spacestep=3);
sol:-plot3d(z,t=0..1000);
sol:-animate(t=1000);