See help(LinearAlgebra[Eigenvectors]): the output of Eigenvectors is a sequence val, vec wher val is a vector containing the eigenvalues and vec the matrix of eigenvectors.
To simplify the output of Eigenvectors do this

simplify([Eigenvectors(L)])

# or

val, vec := Eigenvectors(L):
simplify(vec)

Here is a more synthetic expression of the eigenvalues and eigenvectors:
synthetic.mw

Be aware of undesired effects produced by break, look to this example

restart
for n from 1 to 4 do
  for i from 1 to 2 do 
    if n < 2 then 
      printf("   Within if structure %d  %d\n", n, i)
    else
      break
    end if; 
    printf("Exit if structure %d  %d\n\n", n, i)
  end do;
end do;

   Within if structure 1  1
Exit if structure 1  1

   Within if structure 1  2
Exit if structure 1  2

So just replace

        if Norm(err) <= tol then 
            h := 0.1 * h * (c + 1) * (tol/Norm(err))^(0.2);
        else 
            break
        end if;

by

        if Norm(err) <= tol then 
            h := 0.1 * h * (c + 1) * (tol/Norm(err))^(0.2);
        end if;

Two comments:

Are you sure Digits:=30 is necessary and that  Digits:=10 (default value) would npt ne enough?

Given the large number of points I would prefer one of these two graphics (but it is personal opinion):

numerical_plot_y1 := plot(time_t, numerical_array_y1, style = point, symbol=point, color = blue, legend = ["TFIBF"]):
exact_plot_y1 := plot(time_t, exact_array_y1, style = point, symbol=point, color = red, legend = ["EXACT"]):
display({numerical_plot_y1, exact_plot_y1}, size=[800, 400]);

# or

numerical_plot_y1 := plot(time_t, numerical_array_y1, style = line, color = blue, legend = ["TFIBF"]):
exact_plot_y1 := plot(time_t[1..-2], exact_array_y1[1..-2], style = line, color = red, legend = ["EXACT"]):
display({numerical_plot_y1, exact_plot_y1}, size=[800, 400]);

In Maple 2015 (the version I use) printlevel is an argument of kernelopts, not of interface.
Did you check that?
Note that you can also write printlevel := n directly.

Just a question: Why do you want to change the value of printlevel within your procedure?
If you want to display intermediate results I recommend you to use trace instead:

restart
fac := proc(n::integer)
    if n = 0 then
        1
    else
        n * thisproc(n - 1);  
    end if;
end proc;
proc(n::integer)  ...  end;
fac(5)
                              120
trace(fac):
fac(5)

{--> enter fac, args = 5
{--> enter fac, args = 4
{--> enter fac, args = 3
{--> enter fac, args = 2
{--> enter fac, args = 1
{--> enter fac, args = 0
                               1
<-- exit fac (now in fac) = 1}
                               1
<-- exit fac (now in fac) = 1}
                               2
<-- exit fac (now in fac) = 2}
                               6
<-- exit fac (now in fac) = 6}
                               24
<-- exit fac (now in fac) = 24}
                              120
<-- exit fac (now at top level) = 120}
                              120

Next proceed as described in the attached file:
Ode_New_TWO_phase_mmcdara.mw

Advice...
Odes 2 and 4 both form a subsystem with unknown functions h(Theta) and H(Teta) that you should solve apart.
Note that it can be solved numerically or formally.
Once done, plug this soution into odes 1 and 3 (it's better to have solved the first substystem formally, but you can do the same using the option known [see help(dsolve[numeric]) ] if you solved the first subsystem numerically).

Errors to fix: 

• The first subsystem requires 3 ic/bc and you give 4.
• The first subsystem requires 5 ic/bc and you give 4.

@Prakash J 

Look at this file Difference_mmcdara.mw

@C_R Feel free to delete or reappropriate this comment if you consider that I have interfered inappropriately in your exchanges (I am no longer interested in this race for badges which seems to me to pollute Mapleprimes :-) )

To moderators: this post is deliberately a comment and not an answer, so please don't convert it into an answer. Thanks. 

First point: do not use D[2] (see red highlighted comment in the attached file).

Once corrected, your expression eq is not identically null (for this to happen the constant A, A2, B1, B2, lambda1 and lambda2 must verify some conditions.
This becomes clear if you consider that, for aany strictlyconfinupus function f,  shift(f(n), n) - f(n) is a forward divided difference approximation of diff(f(z), z) at z=n.

verif_May24_mmcdara.mw

Note to the one who converted this comment into an answer: 
Had I wanted to write a reply instead of a comment, I would have done so. It wasn't a mistake on my part, so thank you for not changing again the status of this post.

Download Proposal.mw

[added by moderator, using Maple 2022.2]
(Thanks to the moderator, whoever he/she is)

Download Proposal_2022.2.mw

A slight improvement where unitary gradients are also displayed : Proposal_with_gradients.mw

The attached file contains a numeric exploration of some properties of Eq and an analytic study about its number of real roots.
rho-analysis_mmcdara.mw

As you repliedme once that you wanted analytic results and not results based on numeric simulation or graphic interpretation, this should give you some ideas about how analyzing your equation.

J := n -> int(cos(theta)*LegendreP(n,theta), theta=0..Pi):
J(1)
                               -2
J(2)
                             -3 Pi
J(3);
                               15   2
                          33 - -- Pi 
                               2     
J(n)
      int(cos(theta) LegendreP(n, theta), theta = 0 .. Pi)

Replace

evalf(2*Int(F[m,mu],t=0..7*B[m,mu],method = _Dexp)

by

evalf(2*Int(F[m,mu],t=0..7*B[m,mu],method = _d01ajc)

See the attached file (a print is added in the loop which failed to see the progression of the computations).

LoopError_mmcdara.mw

BTW: I don't understand what you are trying to achieve with your three last lines.

You have primitive pythagorean triples (PPT) and non primitive ones.
For instance [3, 4, 5] is a PPT from which some other triples can be found:

• through permutations you get 6 triples,
• by changing each number to minus its value (if you accept this) you get 8 triples (thus a total of 48 if you use permutations),
• by multiplying each number by any positive integer > 1 you get a new triple.

The key is to find the PPT which are in some sense the generators for other pythagorean triples (PT).

Here is a cowhich build PPTs and some developments to generate induced PT.
PithagoreanTriples.mw

My analysis:

Line 1 , is the equivalent of 

Distrib := unapply(PDF(Normal(average, sqrt(var)), 1.1*x-0.1), (x, average, var))

Line 4  creates a sample (size n) of a random variable with a Truncated Normal Distribution (mean=1, variance=0.399, support=0.8..3).
Statistics doesn't have this kind of distribution but it is very easy to create one.
Sampling it is not a problem neither.
Illustration (code written in a hurry): TruncatedNormal.mw

I don't know what Total means?
Could it be

Mean(func~(RV) /~ Distrib~(RV, 1, 0.399))
# where RV is the name of the sample generated at line (6) of the attached file

If it is so, then the correponding Maple code could be this one:
 

Download TruncatedNormal_2.mw

Here is an example of the uncolding of the prism.
All the steps up to the construction of the Spanning Tree should be generic (at least for a convex polyhedron).
The last one (the drawing of the pattern of this prism) is done by hand.
I didn't try ro make this figure by using the true coordinates of your prism (for the method to apply see Here for instance).

As I don't have time to go any further, I hope that someone else here will be able to take advantage of this code to produce something more substantial.

Download A_first_sketch_for_a_code.mw

Download step_by_step.mw

@acer  and I intertwined: this is basically the same answer but I suggest you to dsolve WITHOUT giving Cl a numeric value.
The solution then depends on t AND Cl.
 

Download Analytic_solution.mw