There's a command that does the image conversion and the file writing in one step: plotsetup. You can use it like this:

currentdir(kernelopts(homedir)); #Set directory where file will go.

ExportJPGMyPlot:= proc(MyPlot, Name::string)
     plotsetup(jpeg, plotoutput= cat(Name, ".jpg"));
     Threads:-Sleep(2); #2 secs to synchronize I/O.
     print(MyPlot);
     plotsetup(default);
     Threads:-Sleep(2);
     MyPlot
end proc:

ExportJPGMyPlot(plot(x^2, x= -2..2), "XSqr");

Note that the API (call sequence, arguments, etc.) is completely different from any other Maple command. It works entirely by side effects which will redirect plots that otherwise would've been displayed in the GUI.

The list of supported image file types is documented (minimally) at ?plot,device. The list of image file types that are available on your system is returned by plotsetup(help). Unfortunately, this latter list is much longer than that documented by the former.

The Threads:-Sleep(2) isn't necessarily needed, and the 2 can possibly be reduced; it's system dependent. I think that the file I/O in this situation isn't synchronized with the Maple kernel. Maybe someone else can confirm or provide further clarification of that.

The command for filtering a set is select. The command fsolve will, by default, return all real roots (and only real roots) for a polynomial. So what you want can be done by

select(x-> g(x,a) > 0, {fsolve(f(x,a))});

To display a 2-D Array A as if it were a Matrix, use <A>.

e:= g^((2*(-sigma+k+1))/(-1+sigma))-tau^2:
ex:= op([1,2], e)/2;
thaw(factor(expand(simplify(e, {ex= freeze(ex)}))));

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)];