The problem is indeed the assignment A1:= A0, which makes the matrices identically the same.

If you want to work entirely in memory, and you have enough memory to do so, replace A1:= A0 with

A1:= copy(A0).

Think of the matrices as boxes. The statement A1:= copy(A0) creates a new box A1 and copies the contents from A0. From that point forward, A1 and A0 can be independently changed. On the other hand, the assignment A1:= A0 says that A1 and A0 are the same box, so any change to the contents of one will affect the other.

If you don't have enough primary memory to store two copies of the matrix, you can store one on disk, like this:

#Create A0.
#Do some calculations with A0.
save A0, "MyA0.m";  #Save A0 for later use.
A1:= A0:  #Uses a trivial amount of extra memory.
#Do some calculations with A1.

save A1, "MyA1.m";
A1:= 'A1': #Optional: Eliminate connection between A1 and A0.
read "MyA0.m"; #Restore A0

If you're interested in real solutions only, then using solve(P, z, parametric, real) gives a little bit more useful information than Kitonum's fsolve. In partcular, it immediately gives a floating-point approximation to the unique real value of w at which the number of real roots changes from 0 to 1 to 2, and it gives RootOf expressions that can be evalfed for those real roots.

I switched the order of your assignments:

restart:
mu:= c^(a-b)/x:
x:= (t^(1-s)+s^(1-2*s))^a:

Now compare mu with eval(mu, 1).

There are many problems:

Why are you looking for complex solutions in the fsolve? If you're expecting nonreal solutions, what do you mean by plotting them?

Are you expecting to get both roots from the fsolve? One is simply the negative of the other. I think that you intend to just use the positive root.

alpha0^2 is so small compared to the other terms in eq1 that your loop that varies alpha0 only makes the roots different after the 40th decimal place.

alpha0 is not a sequence, so alpha0[i] is nonsense. Replace alpha0[i] with 5*(i-1).

Your usage of pi is not an error, but the constant Pi (with a capital P) is already defined, so there's no need for you to set pi:= 3.14.

In the line f5:=, you have a clause Jy= J*sin (theta). You need to take away the extra space between sin and (theta).

Using LinearSolve(A,b) mod 2 won't work: It means to solve the system over the rational numbers and then take that solution mod 2.

The LinearAlgebra:-Modular subpackage is extremely fast. Use it like this:

macro(LA= LinearAlgebra, LAM= LinearAlgebra:-Modular):
vars := [indets(systems, name)[]]:
Ab:= LA:-GenerateMatrix(systems, vars, augmented):
try
     X:= LAM:-LinearSolve(2, LAM:-Mod(2, Ab, integer[]), 1, inplace= false)
catch "matrix is singular":
end try;

The variables can be automatically generated by using the symbol option to Matrix. Then the equations to solve can be automatically generated by converting a Matrix to a list.

b:= Matrix((5,1), symbol= y):
X:= Matrix((2,1), symbol= x):
eval({b,X}, solve(convert(A.b+B.X, list)));

Maple's Optimization package can solve LPs and NLPs. See ?Optimization, ?Optimization,LPSolve, and ?Optimization,NLPSolve.

Setting beta = a+b*I just makes things more complicated. Solve for complex beta.

eq:=
     beta*(1+cos(beta*L)*cosh(beta*L)) +
     I*d*(cosh(beta*L)*sin(beta*L) - sinh(beta*L)*cos(beta*L))
:
L:= 10:  d:= 1:
RootFinding:-Analytic(eq, re= -10..10, im= -10..10);

This'll return 129 roots in under 30 seconds.

If a name x is indexed, then x::indexed will evaluate true, and op(x) will return the sequence of indices. The procedure below is limited to expressions that contain exactly one name.

f:= proc(e::satisfies(e-> nops(indets(e, name))=1))
local x:= indets(e,name)[];
     `if`(x::indexed, {op(x)}, {})
end proc;

e:= (a+b)^3:
[`if`(e::`^`, `$`(op(e)), e)];

The solution to your problem is to use curly braces { } to make a set rather than square brackets [ ] which make a list. In a list, order matters; in a set, it doesn't. The operators union, intersect, minus, and subset will only take sets as operands. (Although it's possible to overload them to work for lists also.) So, make your command

nops({op(x^2*y+x)} ∩ {op(x^2*y+y)});

e:= sqrt(1+diff(f(x),x))*(2+g(x)):
frontend(mtaylor, [e, {diff(f(x),x),g(x)}], [{`+`,`*`,`^`, set}, {}], 2);

My first Reply above was a result of my being thrown off track by your inappropriate use of the word "group", which has a precise mathematical definition. And then I thought that you were looking for a minimal set of generators of the group.

From your followup Reply, and also from VV's clarifying comment, I know that you're looking for a basis (a minimal set of generators) of the subspace of F^(2x2) spanned (or generated) by your given set of matrices, where F is any field containing all of the matrix entries that you present. A basis can be extracted as a one-liner:

map2(Matrix, 2, convert~(LinearAlgebra:-Basis(convert~(S, Vector)), list))^~ %T;

where S is the given set of matrices.

Note that a basis never contains the zero element, which is considered to be generated by the empty set. You'll just need to put that back into the set if you want it there.

Answering the question that is somewhat implied by your title: The set operations for sets with vector (not group) elements that you're looking for are (essentially) Basis, SumBasis, IntersectionBasis, RowSpace, ColumnSpace, and NullSpace (all in LinearAlgebra), but you'll often need to do some fiddling with convert, etc., to get things into the correct format.

Great example. Vote up.

The output's ordering and signage depends on the changing variable's lexicographic order with respect to x0. The y and z come after x0 and all your other changing variables come before x0. The ordering of the the terms in the argument to ln appears to be lexicographic, and the signage (but not the order) of the argument to arctan also depends the variables' lexicographic order. The difference between ln and arctan is, as Tom Leslie notes, that arctan is an odd function. So arctan can do this manipulation. I don't know what purpose it serves to actually do it. The code at showstat(arctan) is fairly brief and fairly easy to understand. You can see some sign manipulations in lines 74 and 75 (Maple 2016) and other places. The result ultimately depends on the what's returned by the built-in sign. Compare sign(x0-x) with sign(x0-y).

You also asked if there was a better way to generate your sequence of variables. Yes:

cat(a..z);