To find diff(f(eta),eta) at the eta-value corresponding to i=3 (as an example) do

i := 'i'; #Necessary, since i has acquired a value.
eval(d1, i = 3);

The other derivatives are similarly found from d2, d3, d4.

Copying your piecewise:

T:=2:
p:=piecewise(t <= 0, 0, 0 < t and t <= T1, dotD21, T+T1 < t and t <= T2, dotD22, t > T2, 0);
p assuming t>T1,t<2+T1;

Question:

Is T1+T <T2? If so the following should be a simpler version of piecewise:
 piecewise(t <= T1, dotD21, t<T+T1,0,t<T2, dotD22, 0);

A semi-concrete version exhibits the problems (I have used the piecewise I suggested):
sol2 := dsolve({eval(de21,{T1=1,T2=4}), de22, ics2});

#####
I suggest using numerical solution:

solnum := dsolve({de21, de22, ics2},numeric,parameters=[T1,T2,U0,L0]);
solnum(parameters=[1,4,1,0.5]);
plots:-odeplot(solnum,[[t,U2(t)],[t,L2(t)]],0..20);

Q:=f(x)*a*b;
#The operands of Q are:
op(Q);
#We see that a*b is not an operand; subs looks for operands (and suboperands):
#So you could do
subs(a=y/b,Q);
#Alternatively, you could use algsubs:
algsubs(a*b=y,Q);

It seems that interpolation fails. Please see the added comment at the bottom.

However, you can save S2 in

restart;
S2 := dsolve([diff(z(x), x$2)+z(x) = -2, z(0)=0, z(1)=1], z(x), 'numeric', 'output' = Array([seq(i/7,i=0..7)])); #You would probably want more values than 8.
plots:-odeplot(S2,[x,z(x)]);
save S2,"G:/temp.m";
restart;
read "G:/temp.m";
S2;
plots:-odeplot(S2,[x,z(x)]);

Also you could if you wish make a procedure that interpolates the values in the saved matrix using CurveFitting:-ArrayInterpolation.

An unrelated comment:
Using unapply to get the procedure for z(x) is like crossing the river to get water:
S1 := dsolve([diff(z(x), x$2)+z(x) = -2, z(0)=0, z(1)=1], z(x), 'numeric', 'output' = listprocedure);
F1:=subs(S1,z(x)); #This is simpler
plots:-odeplot(S1,[x,z(x)],0..1); #Using odeplot is also a good idea.

#####################################################
It seems that somewhere in saving and subsequently reading S1 something happens so that the name of the interpolation procedure `dsolve/numeric/hermite`  no longer refers to a procedure or anything.
But if you do
eval(`dsolve/numeric/hermite`);
in a fresh worksheet, then you see the procedure. You can just select and copy that and paste it into the session after the read statement.

restart;
S1 := dsolve([diff(z(x), x$2)+z(x) = -2, z(0)=0, z(1)=1], z(x), 'numeric', 'output' = listprocedure);
F1:=subs(S1,z(x)); #This is simpler
plots:-odeplot(S1,[x,z(x)],0..1); #Using odeplot is also a good idea.
save S1,"G:/temp.m";

restart;
read "G:/temp.m";
S1;
F1:=subs(S1,z(x));
F1(.12345); #Doesn't work
showstat(F1); #OK
showstat(`dsolve/numeric/hermite`); #Error
##The following procedure is obtained as described above and assigned to `dsolve/numeric/hermite`

`dsolve/numeric/hermite`:=proc (npt::posint, neq::posint, X::Vector, Y::Matrix, YP::Matrix, xv, yv::Vector, L::Matrix, V::array) local n, i, j, m, fac, der, lft, rgt, ipt, lend; ipt := V[1]; lend := V[2]; if lend < 1 then m := `dsolve/numeric/findpoint`(npt, X, xv, V); n := V[1]; V[1] := ipt else m := lend; n := min(ipt, npt+1-m); if n < 1 then error "specified left end point exceeds number of points in mesh" end if end if; lft := evalb(V[3] = true or m = 1 and V[5] = true); rgt := evalb(V[4] = true or m = npt-n+1 and V[6] = true); m := m-1; for i to n do fac := 1; der := 0; for j to n do if i <> j then if j = 1 and lft or j = n and rgt then fac := fac*(xv-X[j+m])/(X[m+i]-X[j+m]); der := der+1/(X[m+i]-X[j+m]) else fac := fac*(xv-X[j+m])^2/(X[m+i]-X[j+m])^2; der := der+2/(X[m+i]-X[j+m]) end if end if end do; if i = 1 and lft or i = n and rgt then L[i, 1] := fac; L[i, 2] := 0 else L[i, 1] := (1-der*(xv-X[m+i]))*fac; L[i, 2] := (xv-X[m+i])*fac end if end do; m := m+1; for j to neq do if lft then yv[j] := L[1, 1]*Y[m, j] else yv[j] := Y[m, j]*L[1, 1]+YP[m, j]*L[1, 2] end if; for i to n-2 do yv[j] := Y[m+i, j]*L[i+1, 1]+YP[m+i, j]*L[i+1, 2]+yv[j] end do; if rgt then yv[j] := Y[m+n-1, j]*L[n, 1]+yv[j] else yv[j] := Y[m+n-1, j]*L[n, 1]+YP[m+n-1, j]*L[n, 2]+yv[j] end if end do; return  end proc;

#Now F1 works:
F1(0.12345);
plot(F1,0..1);

Yesterday I gave you an answer in the link:

http://mapleprimes.com/questions/203185-Implicit-Crank-Nicolson-Finite-Different-Method

My answer (reply) at the bottom was entitled "Latest worksheet".

You haven't seen it perhaps?  (See my comment about email notifications at the bottom).

Actually plots come out of that.

Incidentally, the huge values may be due to wrong signs in the exponential functions (the signs in front of u and s). I think you should look very carefully at what you have copied from wherever it is.

########
Comment about MaplePrimes:

I still find it very strange that even the person asking a question has to subscribe to email notifications. Presumably that person would like to hear about any reply to the question, and if the question receives a flood of replies he could just unsubscribe if he so wishes.

Do you have any reason to believe that a solution exists?

I found a couple of problems with the code you pasted above and have a few comments.

1. after the assignment to n1 there is a line with (presumably) output. I just removed it.

2. In computing the integral for n2 I replaced a by a1 in the range, so I had t = a1 .. a1+tau2-tau1.
After that I did n2:=eval(n2,a1=a);

3. In the code for g1 there appears n1(gam0,gam1) and n2(gam0,gam1). Since n1 and n2 are not defined as functions these should just be n1 and n2.

4. I replaced the occurrences of int in g1 and g2 by Int since I gave up on getting int to work, but would use numerical integration later when values of gam0 and gam1 are supplied as in
evalf(eval([g1,g2],{gam0=0.01,gam1=3.62}));

5. I had no success using fsolve as in
fsolve({g1=0,g2=0}, {gam0=0..230, gam1=0.1..10});

This is not an answer to your question, but I see a few problems in your worksheet.

You consider the interval eta=0..N, where N=10. You have a stepsize of h=0.01. Thus if I understand what you are doing, you should have 1000 steps to take. Yet you seem to be taking only 10.

Secondly, after the loop defining the equations eq[1],..., eq[9] you should be aware that now the loop parameter (i) has a value, namely 10.
That is a problem when you proceed to define eq[10] and the rest, because d1,d2,d3,d4,d12,d14 now has i=10, which you clearly didn't intend.
So right after the loop do this:
i:='i': #Clearing the assignment to i done in the loop

Thirdly, you should think about the meaning of
eq[10] := eval(f[1] = d12) = rhs(bc1[2]);
#Did you not mean
eq[10] := eval(d12,i=0) = rhs(bc1[2]); #or something like it?
#Similarly with eq[11] and eq[12].

#####################
Surely you are aware that the bvp problem is easily solved by dsolve:
restart;
ode1 := (diff(f(eta), eta))^2-f(eta)*(diff(f(eta), eta, eta)) = diff(f(eta), eta, eta, eta)+k_1*(2*(diff(f(eta), eta))*(diff(f(eta), eta, eta, eta))+(diff(f(eta), eta, eta))^2-f(eta)*(diff(f(eta), eta, eta, eta, eta)));
bc1 := f(0) = 0, D(f)(0) = 1, D(f)(N) = 0, (D@@2)(f)(N) = 0;
k_1 := 0.09; N := 10;
res:=dsolve({ode1,bc1},numeric,method=bvp[midrich]);
plots:-odeplot(res,[eta,f(eta)]);

It is often the case that giving an approximate solution helps. Even if that approximate solution is not all that good. If at a first attempt some approximate solution doesn't work, then try another. Remember that Newton's method for systems is used, so starting values are important.
This was my first attempt with an approximate solution for this system and it happened to work:

R := dsolve(eval({Eq1, Eq2, Eq3, Eq4, Eq5, Eq6, bcs1}, B = L[1]), numeric, output = listprocedure,
    approxsoln=[f(eta)=r,F(eta)=0,G(eta)=-r,H(eta)=k,theta(eta)=1,theta1(eta)=0]);

As you will notice the approximate solution is a list of constant solutions, not very sophisticated, but I was inspired by the boundary conditions.

I did this in Maple 18.02. Let me know if it works in Maple 12. I haven't tried it there.

When you look at the graphs it is surprising (maybe) that three of the functions appear to be constant, or almost constant: F, G, and H.

odeplot(R,[[eta,f(eta)],[eta,F(eta)],[eta,G(eta)],[eta, H(eta)],[eta,theta(eta)],[eta,theta1(eta)]],thickness=3,legend=[f,F,G,H,theta,theta1]);
odeplot(R,[[eta,f(eta)],[eta,theta(eta)],[eta,theta1(eta)]],thickness=3,legend=[f,theta,theta1]);
FGH:=[F(eta),G(eta),H(eta)]:
display(Array([seq(odeplot(R,[eta,fcn],thickness=3),fcn=FGH)]));

First remove all output via the menu item Edit/Remove Output/From Worksheet.
Next:
You are obviously not using 1D-math input. If you were there would be no problem: Just copy and paste. Easy.

Although I never use document mode it seems that the following should work:
Select what you want to put in MaplePrimes. Right click. Near the the bottom of the menu you see "2D-Math".
From that chose "Convert To", "1D- Math Input".
Now copy and paste the resulting lines into MaplePrimes.
Don't save your worksheet as the changes will remain in that case.

As an alternative: Upload a worksheet via the thick green arrow in the editor.

To avoid the initial operations necessary for the separation you could just use pdsolve:

restart;
eq1:= diff(T(r,z),r$2)+1/r*diff(T(r,z),r)+diff(T(r,z),z$2);
res:=pdsolve(eq1,HINT=Z(z)*R(r));
eqs:=op([2,1],res);
# Perhaps followed by
eq21,eq22:=op(op~([selectremove(has,eqs,Z)]));

Max:=proc() local num,nam,m, i;
  num,nam:=selectremove(x->type(evalf(x),numeric),[args]);
  if nam=[] then
    m:=-infinity;
    for i in num do
     if i>m then m:=i end if;
    end do;
    m
   elif num<>[] then
    'procname'(procname(op(num)),op(nam))
   else
    'procname'(op(nam))
  end if
end proc;

I had success with other conditions:

restart;
sys:=[diff(C(t,x),t)=In(t,x)*s(t,x)-C(t,x),diff(s(t,x),t)=-In(t,x)*s(t,x),diff(In(t,x),t)=C(t,x)-In(t,x)+diff(In(t,x),x,x)];
ibcs:=[In(t,0)=1,In(t,15)=0,s(0,x)=0.1,C(0,x)=0,In(0,x)=0.9];
pds:=pdsolve(sys,ibcs,numeric,spacestep=0.1);
pds:-value(C)(0.1,.5);
pds:-plot([C,s,In],t=1);
pds:-animate([C,s,In],t=1);

I had to type the equations myself. Not many people (and usually not me either) want to do that. If you paste code as text or upload a worksheet you are more likely to get a response.

1. Notice that if u is an eigenvector so is k*u for k<>0. Thus the results in the two instances can both be correct.
Using v for S and v1 for S1, you could check this for the first column like this
k1:=v[1,1]/v1[1,1]; #The factor
k1*v1[..,1]-v[..,1];
simplify(fnormal~(%));

2. Notice that S is defined before beta has any value assigned to it, whereas S1 is defined afterwards.

Try the following and you may be surprised:

op(S);
%;
op(S1);

3. Be aware that eigenvalues (and thus eigenvectors) may not come in the same order between executions. Although this was certainly true some Maple versions ago, the order may be consistent now,, but I would be cautious and check the eigenvalues for a change in order from last time, or as here between two versions of the same matrix S. The eigenvalues do seem to come in the same order for S and S1 in your case.

4. Different algorithms are most likely used. For the matrix S1 a purely numerical algorithm seems to be used indicated by the fact the 2-norm of all colums is 1.
If you use v for S and v1 for S1 then try
seq(Norm(v[..,i],2),i=1..4);
seq(Norm(v1[..,i],2),i=1..4); #All equal to 1

The exact calculation of the integral 1/pp2 clearly gives you an incorrect answer.

Since pp2 > 0 on the interval considered -int(1/pp2,x=y..7*Pi/16) is clearly negative.
You could try using numerical integration instead:

Tp:=proc(y) if not y::numeric then return 'procname(y)' end if;
-evalf(Int(1/pp2, x = y..7*Pi/16))
end proc;
plot([Tp(y),y,y=0..7*Pi/16],color=blue);

If you are getting anything but {0} then you have a problem.
It could be a simplification problem, i.e. simplify could not see that the expressions you are getting were actually identically zero. To test that you could try inserting (via subs) concrete values for x and for any parameters involved.
Then use evalf. Remember there will be roundoff errors.
Maybe the solution is only valid on an interval about the initial point (if initial values are given).

Simple example:
restart;
sys:={diff(x(t),t)=x(t), diff(y(t),t)=2*x(t)};
sol:=dsolve(sys);
odetest(sol,sys);
#Basically what odetest is doing is this:
eval((lhs-rhs)~(sys),sol);
#However, it can handle all kinds of situations so the procedure is quite lengthy (233 lines):
showstat(odetest);

It could be interesting to see your equations if the answer from dsolve is actually wrong!