It's easy. Let A be the adjacency matrix; so, in this case, A is Y evaluated at the solution, just as you've shown in the PDF. Then

GT:= GraphTheory:
G:= GT:-Graph(A);
GT:-DrawNetwork(G);

The relevant indefinite sums (i.e., for an arbitrary number of terms) for this problem can be done symbolically, and this provides an interesting alternative solution. This type of symbolic analysis is often useful for time series, where the ordinate data are usually simple sequences of indefinite length.

Also, note that there are slightly different formulas for the standard deviation of a sample and of a population. The other Answers and Replys have used the sample formula. It's not clear to me that that's appropriate in this case. So, my code below uses the slightly simpler population formula.

SymbMean:= (f::algebraic, k::(name= range(algebraic)))-> 
   simplify(sum(f,k)/(rhs(rhs(k))-lhs(rhs(k))+1)):
SymbVar:= (f::algebraic, k::(name= range(algebraic)))-> 
   simplify(SymbMean(f^2,k) - SymbMean(f,k)^2):
SymbSD:= (f::algebraic, k::(name= range(algebraic)))->
   sqrt(SymbVar(f,k)):

Usage:
(mu,sigma)=~ Symb||(Mean,SD)(k, k= 1..n); eval(%, n= 20);

                                              (1/2)
                 1     1          1 /   2    \     
            mu = - n + -, sigma = - \3 n  - 3/     
                 2     2          6                
                       21          1     (1/2)
                  mu = --, sigma = - 1197     
                       2           6          
(mu,sigma)=~ Symb||(Mean,SD)(k^3, k= 1..n); eval(%, n= 30);

       1        2            1  
  mu = - (n + 1)  n, sigma = -- 
       4                     84 

                                                         (1/2)
    /     6         5        4         3         2      \     
    \567 n  + 1764 n  + 882 n  - 1764 n  - 1617 n  + 168/     
                 14415          1              (1/2)
            mu = -----, sigma = -- 456873536868     
                   2            84                  

Your instructions say to use a for loop. In Maple, using a for loop for something like this is not necessary, and even a bit unusual and inefficient (unless one is using compiled code or evalhf). But here it is:

Mean:= proc(X::seq(algebraic))
local S:= 0, x;
   for x in X do S:= S+x od;
   S/nargs
end proc
:
Mean(seq(k^3, k= 1..30));

Since the variance and standard deviation can be expressed as simple functions of the mean, using a for loop for those is just ridiculous.

Here's a little procedure that creates a new random variable from any real-valued algebraic combination of any number of any independent random variables with finite support (hence discrete). The new random variable will be treated with the same respect and efficiency by the Statistics package as the originals. The underlying random variables need not be identically distributed, nor equiprobable, nor have their support be a range of integers. So, if I've read the other Answers and Replys correctly, I believe that this does everything that they do and more.

restart
:
St:= Statistics
:
CombineDiscrete:= proc(R::list(RandomVariable), f)
uses S= Statistics, It= Iterator;
local r, X, t, P:= table(sparse), nR:= nops(R), fX;
   for X in 
      It:-CartesianProduct(
         seq([solve~(op~(indets(S:-PDF(r,t), specfunc(Dirac))), t)[]], r= R)
      )
   do
      P[(fX:= f(seq(X)))]:= P[fX] + mul(S:-Probability(R[r]=X[r]), r= 1..nR)
   od;
   P:= [entries(P, 'pairs')];
   S:-RandomVariable(EmpiricalDistribution(lhs~(P), 'probabilities'= rhs~(P)))
end proc
:  
Dice:= 'St:-RandomVariable(DiscreteUniform(1,6))' $ 2:
Game:= CombineDiscrete([Dice], (x,y)-> piecewise(x=1, -4, y=1, -4, abs(x-y)) ): 
St:-Mean(Game);
                               -1
                               --
                               9 
St:-Variance(Game);
                              620
                              ---
                              81 
St:-Probability(Game=3);
                               1
                               -
                               9
St:-Probability(Game > 2);
                               1
                               -
                               6
seq(St:-Sample(Game, 32));
  1., 2., 1., 3., -4., 2., -4., 2., 2., -4., 0., 1., 0., -4., 1., 
  1., 2., -4., 2., 1., 3., 1., -4., 3., 2., 2., 0., -4., 3., 1., 
  2., 4.
S:= CodeTools:-Usage(St:-Sample(Game, 10^6)):
memory used=15.27MiB, alloc change=15.27MiB, cpu time=47.00ms, real time=55.00ms, gc time=0ns

Doing symbolic computation by literal string manipulation is tough. Most things can be and should be done by checking the type of an expression. If a polynomial has been sucessfully factored, then its type will be `*` or `^` and not `+`. So, the if in your inner loop can be

if not temp4::`+` then found one fi;

It is possible to solve this problem in a manner akin to Rouben's Answer without knowing anything about Fibonacci numbers, without doing any mathematical induction, without knowing which position in the matrix is largest, and without putting an artificial upper limit (such as 30) on n; and in a way that can be generalized to similar matrix problems. I was intending to post this even before seeing Rouben's Answer anyway. Like this:

restart:
A:= <1, 1; 1, 0>:  M:= Matrix((2,2), symbol= a):
S:= rsolve(
   {
      seq(apply~(M,n+1) =~ A.apply~(M,n)), #recurrence equations
      seq(apply~(M,1) =~ A) #initial conditions
   }, 
   {seq(apply~(M,n))} #functions to solve for
);
N:= min(seq(fsolve(eval(abs(v(n)), S)=2019, n= 1..infinity), v= seq(M)));
                        N := 16.48726057
is(2019::RealRange(max(A^floor(N)), max(A^ceil(N))));
ceil(N);
                         17

In a similar solution technique, Maple can compute the general entry of A^n (for symbolic n) using eigenvalue methods (the specific command that uses the eigenvalue methods is LinearAlgebra:-MatrixFunction):

restart:
A:= <1, 1; 1, 0>:
min(fsolve~(abs~(LinearAlgebra:-MatrixFunction(A, x^n, x)) =~ 2019, n= 1..infinity));
                   16.48726057 

Now here's a completely different solution technique. The above methods involve no guesswork. But if we're willing to make a very small and reasonable guess---that max(A^n) is approximately an exponential function of n---then a very efficient solution ensues. We compute A^(2^n) by repeated squaring. Once the target value is exceeded, we interpolate an exponential function between the last two computed powers.

MaxMatEntry:= proc(A::Matrix, M::positive)
local `A^(2^n)`:= copy(A), k, oldmax, newmax:= max(A), a, b, n;
   for k do
      oldmax:= newmax 
   until (newmax:= max((`A^(2^n)`:= `A^(2^n)`^2))) > M;
   eval(n, fsolve({a*b^(2^(k-1))=oldmax, a*b^(2^k)=newmax, a*b^n=M}))
end proc
:
MaxMatEntry(A, 2019);
                          16.48726041

Note that this agrees with the previous answers to 6 decimal places, even though we only care about the integer part.

The command max will return the maximum entry of a Matrix (or any other container (set, list, Vector, etc.)).

Although it won't make a measurable difference for this small problem, here's a trick that's lets you avoid recomputing the powers. We use the obvious identity A^n = A^(n-1) . A:

restart:
A:= <
   1, 1;
   1, 0
>:
`A^n`:= copy(A):
for n from 2 do until max((`A^n`:= `A^n`.A)) > 2019:
n;

Notes:

By using the "name quotes" ``, a variable's name can be made from any expression, and we can choose these in way meaningful to a reader. So, `A^n` is just a variable name.

The copy is necessary because we need a new Matrix to be created. B:= A just creates a new pointer to A rather than a new container. 

I think that you're looking for the Cartesian product A x B x C. Here are two ways to get that in Maple:

A:= [1,5,7]:
B:= [23,56,12]:
C:= [1,3]:
CartProd1:= (L::seq({Vector, list}))-> 
   [seq(Array(`..`~(1, numelems~([L]))[], ()-> `?[]`~([L], `[]`~([args]))))]
:
CartProd1(A,B,C);

CartProd2:= (L::seq({Vector, list}))->
   [seq([seq(p)], p= Iterator:-CartesianProduct(L))]
:
CartProd2(A,B,C);

The two methods give their results in different orders.

It doesn't make sense to use a set for this because you can't control the order of the elements of a set. But you can use a list or a Vector. 

BoundedSum:= proc(L::{list, Vector}, U::realcons)
local k, S:= 0;
   for k to numelems(L) do until (S:= S+L[k]) > U;
   L[..k-1]
end proc:

A:= [20,15,10,30,46,78]: #Use [ ] or < >, not { }.

BoundedSum(A, 50);

The above works in Maple 2018. If needed, I can easily retrofit it for earlier Maple.

I cannot fault Tom Leslie for having given up trying to understand what you were trying to do; your number of errors is extreme. However, having read through your code numerous times over the past 17 hours, I figured out that you are trying to find a degree-12 power series solution to the IVP

   y'''(x) = -1/2*y(x)*y''(x), y(0) = 0, y'(0) = 1, y''(0) = A.

That can be done like this:

Order:= 13:  # "environment" variable for degree of "order" term of power series 
S:= dsolve(
   {diff(y(x),x$3) = -1/2*y(x)*diff(y(x),x$2), y(0)=0, D(y)(0)=1, (D@@2)(y)(0)=A},
   y(x), 'series'
);
G:= unapply(convert(S, 'polynom'), x); 

Pay close attention to the where I have placed the parentheses! Also, parentheses and square brackets are not interchangable. On the other hand, the spacing and line breaks are a matter of personal style and can be changed however you desire.

Depending on the size of B and the uniqueness of its elements, there can be a considerable benefit to converting B to a set first:

remove(member, A, {seq(B)})

Let n be the number of elements of A, let m be the number of elements of B, and let q be the number of distinct elements of B. Then (from theoretical considerations alone) the number of operations to perform remove(member, A, B) is proportional to n*m, and those to perform remove(member, A, {seq(B)}) is proportional to n*log(q)+m*log(m).

The lowercase e can only be used in explicitly specified floating-point constants. If you need to use a variable, then use 10. as a base, as shown by Acer.

The L^p norm of a function f (x) over an interval a..b is an integral---int(abs(f(x))^p, x= a..b)^(1/p)---for 1 <= p < infinity. The limit of this as p--> infinity is the "infinity norm", which (for continuous f) is also the maximum of the absolute value. So, Tom Leslie's Reply is incorrect about the L^2 norm that you asked about. He's just inadvertently recomputed the infnity norm (which can be clearly seen from his numeric ranswer).

Here's a very short procedure to compute it for any p, including infinity:

LpNorm:= (f::algebraic, xab::(name= range(complexcons)), p::positive)->
   `if`(p=infinity, numapprox:-infnorm(f, xab), int(abs(f)^p, xab, numeric)^(1/p))
:

Using it for your example:

Digits:= 15: #Good enough for most purposes.
exact:= x^4-x^3:
app:= -2.220446049250313e-16+.19676728660169784*x+4.6491983150099205*x^2
         -3.9010072748158082*x^3+2.3619056303875987*x^4-.33198038154796344*x^5:
map[3](LpNorm, exact-app, x= 0..2, [1, 2, infinity]); 
      [6.36243961193564, 5.55700549477000, 6.94938751138337]
 

You need to also save the whole solution (i.e., the direct output of your dsolve(..., numeric) command). If you save a snippet of code that refers to some other variable that you've defined, then you need to save that other variable also. Since X_br refers to solution, that is the case here. Note that you can save any number of variables together in a single file using a single save command. Here's how to do it:

restart
:
solution:= dsolve({diff(xbr(t),t$2) = -xbr(t), xbr(0)= 0, D(xbr)(0)= 1}, numeric):
X_br:= tau-> eval(xbr(':-t'), solution(tau)):
save X_br, solution, "/Users/carlj/desktop/NumSol.mpl"
:
#Now move to your plotting worksheet.
restart
:
read "/Users/carlj/desktop/NumSol.mpl":
plot(X_br, 0..7);

There are three important, significant differences between my code and yours:

1.  I do not use output= listprocedure ^[*1].
2. My X_br is a procedure; yours is an expression.
3. I include solution in the save command.

And there's one minor significant difference: I used the procedure-form plot command

plot(X_br, 0..7);

rather than the expression-form

plot(X_br(t), t= 0..7);

If you have some good reason for needing the expression form, that's still possible. Let me know.

The only difference that using a ".m" internal-format file makes is that Maple is a tiny bit more efficient working with those files. Sorry about my incorrect previous Answer regarding ".m" files.

^[*1]Actually, I never use output= listprocedure. I think that it's a shoddy programming practice, and there's actually nothing that's made easier by it.

You are right to be suspicious: The way that you propose to iteratively modify sets is quite inefficient. What you want to do can be done efficiently^[*1] by less than a quarter of a line of code:

{Ug3[], seq(seq(u..6, -6), u= Ug3)};

Note that this doesn't require the unwanted numbers to be removed because they aren't generated in the first place.

You asked if sets are efficient. Sets and lists are very efficient if you don't try to repeatedly modify them. If you need to repeatedly or iteratively modify the structure, use a Vector or table, and if you ultimately need a set or list, convert it to such when you're done with the modifications.

^[*1]Tom Leslie's one-liner is slightly more efficient than my code above.

Make the filename "filename.m" instead of "filename.mpl". The .m form is called "internal format". In addition to the code of the procedure, this format also puts in the file all the global variables needed to  execute the procedure correctly. 

[Edit: I was wrong about this. Ignore this Answer.]