You forgot colon in assigning L. Also evalb command can be omitted:

f:={1,2,3,4}:

h:={1,2,4,5}:

L:=seq(i,i=1..4):

for i in L do

if  f[i]=h[i]  then

print(f[i]);

end if;

end do;

One way, you can do some calculations in for loop and in each step assign a name to your result.

An example:

for a from 0 to 3 by 0.5 do  # For each value of parameter a we find the list of roots of the equation  x^3-a*x+1=0

R[a]:=[fsolve(x^3-a*x+1)];

od;

Another way to use  seq  command (a solution of the same problem):

R:=seq([fsolve(x^3-a*x+1)], a=0..3, 0.5);

Reference to these results you have in different ways: in the first case, the indices coincide with the values of the parameter a, in the second case as the index, you just specify the number of the solution.

Edited.

Your points do not lie on one straight line. I guess that you want to find a straight line using the least squares method. Here is the solution together with visualization:

restart;

LogA[t] := [- 0.097,  -0.20, -0.22, -0.25, -0.30 ,-0.40, -0.45, -1.01]:

T := [0, 20, 40, 60, 80, 100, 120, 140]:

A:=plot(T, LogA[t], style=point, color=blue, symbolsize=17): 

Eq:=CurveFitting[LeastSquares](T, LogA[t], t);   # Equation of straight line

B:=plot(Eq, t=T[1]..T[-1], color=red, thickness=2):

plots[display](A, B);

I never use verify  command; expand command solves the problem:

expand(cos(u)+sin(u) - sqrt(2)*cos(u-Pi/4));

                                           0

In total we have 16 solutions. Displayed 8 solutions that all numbers are positive:

sols := [isolve(sum(2*(m+i-1)+1,i=1..n) = 3375)]:

S:=select(hastype, sols, identical(n) = posint):

nops(S);

for s in S do

if eval(m,s)>=0 then print([seq(eval(2*m+1+2*(n-1),s[1]), n=1..eval(n,s[2]))]) fi;

od;

Without loss of generality, we can assume that  L>0 :

x1:=(sin(beta*x)*cos(m*Pi*x/L))/(1+alpha*x);

int(x1, x=0..L) assuming  L>0, alpha>0; 

If  alpha<0 , then in the integration interval can be a point of discontinuity of  x1 . This case requires separate consideration.

Your problem can be solved by  LinearAlgebra[Eigenvectors]  command. Below we see that  w1=w  and  Q11=Q

Download 1234567.mw

A corrected version (in several places the multiplication sign was missed or was written incorrectly):

Download 1234.mw

First way:

b:=y=x^3;

op(remove(has, [solve(b, x)], I));

Addition.

A more reliable way (the example of Preben):

b:=solve(x^3-3*x+1=0, x):

map(t->Re(t), remove(t->is(Im(t)<>0), map(evalc,[b])));

To get a predicted response in the presence of the parameters (in the equation  y=x^3  variable  y  is a parameter) necessary to impose a restriction on parameters:

b:=solve(y=x^3, x):

assume(y>0):

map(t->Re(t), remove(t->is(Im(t)<>0), map(evalc,[b])))[];

                                       y^(1/3)

or

b:=solve(y=x^3, x):

assume(y<0):

map(t->Re(t), remove(t->is(Im(t)<>0), map(evalc,[b])))[];

                                    -(-y)^(1/3)

Edited.

Maple doesn't solve such equations in matrix notations. If the matrices  B  and  C  are known, the equation can be easily reduced to a system of polynomial equations that Maple can solve numerically.

An example:

A:=Matrix(2, symbol=x);

B:=Matrix(2, [-1,2,3,4]):

C:=Matrix(2, [-1,0,4,-5]):

Sys:={seq(seq((A^2+B.A+C)[i,j], j=1..2), i=1..2)};

evalf(allvalues(solve(Sys)));

The equation has no solution. 

Here are my arguments:

1. Your equation does not contain the unknown function  y(x) (only its derivatives) and your initial conditions  do not contain any conditions on  y(x)  itself. Therefore, if the equation has a solution  y(x) , then any function  y(x)+constant would also be a solution. 

2. I added the condition  y(0)=0  and dropped the condition  D@@2)(y)(0) = 0 .  We got  a unique solution. But the plot of the second derivative clearly shows that  (D@@2)(y)(0) <> 0 :

eq := diff(y(x), `$`(x, 3))-(diff(y(x), x))*(diff(y(x), x))+1 = 0;

cond := y(0)=0, (D(y))(0) = 0, (D(y))(1) = 1;

de := dsolve({cond, eq}, y(x), numeric);

plots[odeplot](de, [x, diff(y(x),x,x)], x=0..1);

Addition. Without any plot:

de(0):

                      [x = 0., y(x) = 0., diff(y(x), x) = 0., diff(y(x), x, x) = 1.39362859931605]

Maybe Lyapunov function will help you.

An example for lists:

M:=[seq(i^2, i=1..20)]:  # list 1

N:=[seq(i^3, i=1..20)]:  # list 2

plot(M,N);  # or

plot(M,N, style=point);

The same example for vectors:

M:=Vector(20, x->x^2):

N:=Vector(20, x->x^3):

plot(M,N);  # or

plot(M,N, style=point);

All the types of arrows you can see in the help on  DEtools[DEplot].

For example:

DETools[phaseportrait]([D(x)(t)=-0.4*x(t)+(0.5+4*x(t))*y(t),D(y)(t)=0.4*x(t)-(4.5+4*x(t))*y(t)],[x(t),y(t)], t=0..100, [[x(0)=1,y(0)=0.0]], x=0..1,y=0..0.1, stepsize=0.01, scaling=UNCONSTRAINED, linecolour=BLACK, dirgrid=[17,17], linestyle=1, arrows=mediumfill, axes=BOXED);

If we assume that the parameters  r  and  theta  are real numbers,  evalc  command calculates the initial expression to the complex number written in algebraic form:

sqrt(1-r^2/exp(2*I/theta));

simplify(evalc(%));

 Edited  (correction):

simplify(evalc(sqrt(1-r^2*exp(2*I*theta))));