X3 := Groebner:-LeadingMonomial(g3, plex(x,y,z,w)); 

However the leading term in lex order with x > y > z > w is -x*z not y^2.

I think the coefficients are large and it is taking a very long time to compute them.  You can try to compute the Groebner basis modulo a small prime instead and perhaps that will give you some information.  Set infolevel[Groebner] := 5: to watch the progress:

infolevel[Groebner] := 5:

G := Groebner:-Basis(LL, 'tord', characteristic=65521):

Unfortunately solving polynomial systems tends to produce exactly this kind of blowup, and there is often little anyone can do.  The algorithms generally have to conserve as much memory as possible in order to produce any solution at all

If you post the system we can try to solve it, but it could also be the case that your problem is too hard at present.

This is an area of active research so I encourage you to post the system.

It's possibly a bug or a bad case for `evala/Factors/efactor` which is the final factorization method used by Maple with algebraic extensions.  You can see the code:

kernelopts(opaquemodules=false): interface(verboseproc=3):

eval(`evala/Factors/efactor`:-diophant);

The help page for ?evala,Factors says you can try different factorization methods, but I don't see this in the code.  I think the documentation is out of date.

evala(Factor(f,{I,sqrt(3)}), trager);

evala(Factor(f,{I,sqrt(3)}), lenstra);

You should try to post some source code that we can run.

f := randpoly([x,y,z]);  # expression

fn := "filename.txt": seq(fprintf(fn,"%a\n",i), i=f): fclose(fn):

Also of interest, number of objects and memory usage by type:

kernelopts(memusage);

subs can't isolate A^2+1 in the outermost level.  A quick fix is to substitute for A^2.

f:= A^2-2*sqrt(A^2+1)+1;

subs(A^2=M^2-1, f);

Then you may want to simplify it, etc.

Maple automatically applies the distributive law "a*(x+y) = a*x+a*y" for numbers.  For programming purposes, e.g. fraction free algorithms, you can get a numerator and common denominator from the numer and denom commands.  E.g.:

f := 1/3 + x/5;

n := numer(f);

d := denom(f);

However if you form the fraction n/d it will again produce 1/3+x/5.

It's hard for us to try with just the image you posted.  Can you upload the worksheet?  Also (and it looks like you are doing this) it is much easier to simplfiy(a-b) to zero than it is to simplify(a) and simplify(b) to the same thing.  Another thing you could try is the is() command, as in is(a-b=0);  It returns true (probabilistically true), false (definitely false), or FAIL (doesn't know).

Yes but having a loop inherently exposes shared variables.  You need to make the body of the loop into a procedure where each iteration can run independently, then you can call Threads:-Seq or Threads:-Map.

For the purposes of solving, writing v(t) will define v as a function of t.  E.g.:

diff(f(t)*g(t),t);  # returns f'(t)*g(t) + f(t)*g'(t)

A 10000 x 10000 matrix with hardware floating point numbers is 75GB.  As a bit matrix it would be 1.16 GB.  You might want to look into the Bits package.

gcd and gcdex are for polynomials, igcd and igcdex are for integers.

a := 35:

b := 21:

g := igcdex(a,b,'r','s');

r; s;

g = a*r + b*s;

This is certainly possible using a C string and OpenMaple.  It might help to use a backslash to separate lines of Maple code in the string.

http://www.maplesoft.com/support/help/Maple/view.aspx?path=OpenMaple/C/MapleEval

@NoThik 

Alternatively, you could create the equations as polynomials rather than use matrices.  E.g. instead of creating a matrix equation A*X = B, where A := [[1,-1],[1,1]]; and B := [1,-1]; you can create the polynomials {x1-x2 = 1, x1+x2 = -1}.  Internally, Maple converts to this representation to solve sparse linear systems over the rationals.  You can watch what the solver is doing by setting infolevel[solve] := 5;  Here is Markiyan's example, converted to polynomials:

with(LinearAlgebra):
m,n := 2000,2000:  # rows and columns
A := RandomMatrix(m, n, density = 0.01, generator = rand(-1 .. 1));
B := RandomVector(m);

var := [seq(cat(x,j),j=1..n)]:
sys := [seq( add(A[i,j]*var[j], j=1..n) - B[i], i=1..m)]:
infolevel[solve] := 5:
sol := SolveTools:-Linear(sys,var):

If your system is highly structured this solver is fast.  It will be called if you use LinearAlgebra:-LinearSolve as well, but there will be conversion overhead in both directions.