Is this what you're looking for?

restart:
A:= <
   1.364, 0.64765768e-1  ; 
   2.05,  -.182176113    ; 
   2.664, -0.13914542e-1 ; 
   2.728, 0.2193938e-1   ; 
   4.092, -0.18349139e-1 ; 
   4.1,   -.312968801    ; 
   5.328, -0.1819837e-2  ; 
   5.456, -.28840961     ; 
   6.15,  -.57076866     ; 
   7.992, .175022254
>: 
F:= unapply(CurveFitting:-LeastSquares(A, x), x):
P:= plot(
   [F, A], 0..8, color = [red, blue], style = [line, point],
   legend= ["Метод наименьших квадратов", "Экспериментальные данные"], 
   legendstyle= [font = ["Roman", 15]],
   labels= [typeset(d, ", ", `мм`), typeset(ln(`I__ `/I__0))], labelfont = ["Roman", 15],
   axesfont= ["ROMAN", "ROMAN", 15], linestyle = [solid], symbolsize = 20, 
   title= "Определение линейного коэффициента поглощения", titlefont = [Roman, bold, 20]
);
X:= A[..,1]:  Y:= A[..,2]:  `Y^`:= F~(X):
plots:-display(
   [
      P,
      Statistics:-ErrorPlot(
         `Y^`, xcoords= X, yerrors= abs~(Y-`Y^`), markers= false, thickness= 0
      )
   ],
   axes= normal, gridlines= false
);

Rather than putting an upper bound on k in the for loop, you could put an upper limit on the time your willing to spend on the computation. Here's my version:

Wrapped_prime:= proc(p::prime, {base::posint:= 10, filter::anything:= (()-> true)})
global lastK;
local pow:= base^trunc(ln(p)/ln(base)), b2:= base^2, k, m:= p;
   lastK:= 0;
   for k do
      #Recursive computation of repunit-wrapped number:
      pow:= b2*pow;  m:= pow + base*m + 1;
      if filter(k) then 
         if isprime(m) then return k fi;
         lastK:= k
      fi
   od
end proc
:     

Since repunits exist in any base, not just base 10 (Mersenne numbers being the most-famous example), the procedure will work in any base, defaulting to 10. The keyword parameter filter lets you specify a filtering of the k values, for your example, it's filter= (k-> irem(k,6) = 0). The  global lastK lets you check what was last value of k for which the wrapped number was verified composite by isprime.

The first argument to timelimit is the number of seconds of CPU time to expend before giving up.

timelimit(5, Wrapped_prime(397, filter= (k-> irem(k,6)=0)));
Error, (in isprime) time expired
lastK;
                              1656
timelimit(5, Wrapped_prime(397, base= 5));
                               2
timelimit(5, Wrapped_prime(397, base= 2));
                               24

Your system of equations contains numerous terms of the form

sum(FyjR, j = 1 .. m),

which makes it seem as if you believe that FyjR can be indexed by j; it can't be. That sum evaluates to simply m*FyjR.

Maple's mod is built to do computations in higher algebraic structures (such as rings of polynomials and matrices) built on top of the ring of integers modulo m. Thus, unevaluated variables such as your n are treated as transcendentals (or "dummy variables") in a polynomial, so n mod 2 is truly equal to n. In your case, you want n to represent an integer (which presumably will be supplied later), not a transcendental. In that case you should use irem(n,2) instead of mod(n,2).

While I strongly recommend that you use irem instead mod for this situation, Tom Leslie shows that you can use mod(n,x) if both n and x are unevaluated and integer values will be supplied for both simultaneously, or at least n is supplied before x.

Okay, you're running the command-line version of Maple interactively. Let's suppose that the file that you created with Vim is named "C:/dir1/dir2/MyCode.mpl". Then, at the ">" cursor/prompt, enter

read "C:/dir1/dir2/MyCode.mpl";

The fact that your function involves exp is insignificant. You need numeric values and/or ranges for the variables t, a, and eps or else clearly you'll be "unable to evaluate the function to numeric values". Also, we must have -2 < a < 2 to avoid complex values. 

Try this command:
Explore(
   plot(2*cos(t)/sqrt((1+exp(-eps*t))*(4/a^2-1)), t= -Pi..Pi, view= [DEFAULT, -2..2]), 
   a= 0.1..1.9, eps= 0.0..1.0
);

This'll let you control the values of a and eps with sliders.

Here's what your final plot should look like, more-or-less (the colors, arrow size, etc., are up to you):

This shows that the Euler's method solution is very wrong! That's the point of the problem. The step size is too large.

If you didn't get a plot like above, show your code, and we'll suggest corrections.

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.