The easiest way is to write a recursive program that uses branch and bound. The program should generate the possible values of a variable x and call itself recursively to find values for the remaining variables that sum to n-x. My program is below, the equation x+y+z would be represented by a list of coefficients [1,1,1].

recsolve := proc(L, n)   
  if nops(L)=1 then
    `if`(irem(n,L[1],'q')=0, [args[3..-1],q], NULL)
  else
    seq(procname(L[2..-1], n-i*L[1], args[3..-1], i), i=0..iquo(n,L[1]))
  end if
end proc:

#example:  solve for [1,2,3].X = 5
recsolve([1,2,3],5);

This comes from a program I wrote some time ago that found all positive integer vectors whose dot product with a given positive integer vector was in some non--negative range. It essentially solves a linear diophantine inequality in the case where the feasible region is bounded. You could probably adapt it to solve a system of linear diophantine inequalities, but I never got around to doing that.

recrange := proc(L, a, b) local i, r;
  if nops(L)=1 then
    seq([args[4..-1], i], i=iquo(a,L[1],'r')+`if`(r=0,0,1) .. iquo(b,L[1]))
  else
    seq(procname(L[2..-1], max(a-L[1]*i,0), max(b-L[1]*i,0), args[4..-1], i), i=0..iquo(b,L[1]))
  end if 
end proc:

# example: solve [1,2,3].X for the range 3..5 
recrange([1,2,3], 3, 5);
seq(inner(i,[1,2,3]), i=%); # check

In general you can solve multivariate polynomial systems using resultants, triangular sets, or Groebner bases. The resultant of two polynomials f and g with respect to a variable x eliminates x, and for polynomials in two variables this is really the fastest method. Example:

f := x^2*y - 1; 
g := x*y^2 + y - 1; 
r := Resultant(f,g,x) mod 11; # polynomial in y
Factor(r) mod 11;  # roots in y

You can read off the roots in y and for each root substitute it into {f,g} to obtain polynomials in x. Their gcd contains the corresponding values (if any) of x. In the example above, the resultant has two factors y and y+7 which have integer solutions mod 11. Lets try:

S := subs(y=0, {f,g}) mod 11;  # produces non-zero integer, no solution for that
S := subs(y=-7, {f,g}) mod 11; # produces polynomials in x
S := Gcd(op(S)) mod 11;  # produces x+5

Then {y=-7, x=-5} is a solution mod 11. In fact, it is the only solution using the integers 0 through 10. The remaining factor of the resultant introduces an algebraic extension, which you may not be interested in. It is basically introducing the root of a polynomial in y to solve the system, similar to how sqrt(2) could be introduced to "solve" y^2-2=0 over the rational numbers. Keep in mind that resultants may produce extraneous roots, such as the root y=0 above. You have to check the roots against the original system. If you're just interested in a quick and dirty solver for small problems, try the msolve command. For example:

msolve({f,g},11);  # produces {y=4,x=6}

This command is somewhat limited to small systems and integer roots. If you need to solve larger systems with more variables or more general coefficients let me know and I can describe some of the other methods available in Maple.

The problem in your code is that you don't really have a handle on the variable names you want to unassign. Your code references them indirectly as lhs(x)[i], which can be evaluated to the value of the name, but not very easily to the name itself. When programming in Maple, it's usually best not to think like C. Try the following:

AssignVariables := proc(x) 
  local i, n, v, r; 
   v := lhs(x);
   r := rhs(x);
   n := op(1,v);
   for i to n do
     v[i] := r[i];
   end do;
end proc:

Because x is a name/table/rtable/Matrix/Vector you can assign to it in a procedure. Don't worry about using the assign and unassign commands, they are for situations where you can't use the assignment operator (like in a seq, add, or mul).

Put the code into a text file and edit it with your favorite text editor. Read the file into Maple using the read command.

Maple can not assign to procedure parameters, but there are a few ways of working around this. Note that Maple automatically returns the last value of a procedure, so for simple functions you don't need an explicit return. 1) make the procedure use a global variable instead of a parameter

x := 0:
inc := proc() global x; x := x+1; end proc: # returns x+1
for i from 1 to 10 do inc(); end do;

2) return the modified parameter sequence and assign it

fib := proc(f1,f2) f2, f1+f2 end proc; # returns (f2, f1+f2)
f1,f2 := 1,1;
for i from 1 to 10 do f1,f2 := fib(f1,f2); end do;

3) pass the parameter's variable (unevaluated) and assign to it

x := 0;
inc := proc(x,y) y := x+1; end proc: # x=value, y=variable
for i from 1 to 10 do inc(x,'x'); end do;

4) pass the parameter's variable and evaluate it when reassigning

x := 0;
inc := proc(x) x := eval(x)+1; end proc:
for i from 1 to 10 do inc('x'); end do;

Methods 3) and 4) look a lot like passing a pointer in C, however I implore you not to use them, especially 4). Maple's eval command is not the same as pointer dereferencing, and depending on what x is it can be very expensive. Method 3) can be useful when you need to modify something in a deeply nested seq, however that's not something you are likely to need often. Method 2) is the best in simple situations. It's the easy to understand and easy to debug, and it uses a great feature of the language. Method 1) is probably the most widely used, because x typically isn't global variable. It can be local to some procedure or module, and it will be lexically scoped automatically. Here is an example using a module:

counter := module()
  local x;
  export inc;

   x := 0;  # initialize
   inc := proc() x := x+1; end proc;
end module:

counter:-inc();
counter:-inc();

We can make this module behave like a procedure as well, by renaming inc to ModuleApply. A procedure with this special name is called automatically in the situation below:

counter := module()
  local x, ModuleApply;

   x := 0;  # initialize
   ModuleApply := proc() x := x+1; end proc;
end module:

counter();
counter();

The end result is a black box state machine. It's not really a black box though, you can crack it open with a kernel option:

kernelopts(opaquemodules=false):
eval(counter);
eval(counter:-x);
eval(counter:-ModuleApply);

Maple 10 can solve this in ~15 seconds:

TIMER := time(): 
sys, vars := {a1*alpha+b1 *beta=a1^2*alpha+2*a1*a2*gamma+a2^2*gamma*delta/beta,a2*alpha+b2* beta=a1^2*beta+2*a1*a2*delta+a2^2* (beta*gamma-delta*(alpha-delta))/beta, a1*gamma+b1*delta=a1*b1*alpha+a1*b2*gamma+a2*b1*gamma+a2 *b2*gamma*delta/beta,a2*gamma+b2*delta=a1*b1*beta+a1*b2* delta+a2*b1*delta+a2*b2*(beta*gamma-delta*(alpha-delta))/ beta, a1*gamma*delta/beta+b1*(beta*gamma-delta*(alpha-delta))/ beta=b1^2*alpha+2*b1*b2*gamma+b2^2*gamma*delta/beta,a2* gamma*delta/beta+b2*(beta*gamma-delta*(alpha-delta))/ beta=b1^2*beta+2*b1*b2*delta+b2^2*(beta*gamma-delta*( alpha-delta))/beta},{a1,a2,b1,b2}:

# rewrite as polynomial equations
# replace gamma=0.5772156649 by a harmless name
sys2 := subs(gamma=GAM, map(expand@numer@(lhs-rhs), sys));

# factor polynomial system
with(PolynomialIdeals):
J := PolynomialIdeal(sys2, variables=vars):
P := Simplify(PrimeDecomposition(J)):

# generic solutions
S := seq(solve(Generators(i), vars), i=P);
time()-TIMER;

This method is usually faster than solve when there are a finite number of generic solutions.

The numtheory[cfrac] command will probably do what you want. with(numtheory): p := cfrac(Pi); p1 := cfrac(p);

I requested this feature some time ago. You can't do it in Maple 10.

> map(normal, expr); Where expr is your expression. This will apply the normal command to each term in the sum. You may also want to try: > map(factor, expr); However that will factor (1-e^2) as (1-e)*(1+e).

The command you want is "simplify with side relations", where the second argument is a set of relations. For example:

simplify(-1/4*rho*Pi*a^4+1/4*rho*Pi*b^4,  {-rho*a^2*Pi+b^2*rho*Pi=mass});

See also the help page ?simplify,siderels. Because mass appears on the right hand side of an equation, the expression will be rewritten in terms of mass as much possible. You can also choose the order explicitly with a third argument:

simplify(-1/4*rho*Pi*a^4+1/4*rho*Pi*b^4,  {-rho*a^2*Pi+b^2*rho*Pi=mass}, [rho, a, b, mass]);

This will eliminate rho as much as possible, followed by a, then b. If you want to make one of the variables invertible (ie: it can appear in a denominator), leave it out of the list.

Setting infolevel[solve] := 5; will output information about the solving process, but unless the input is a linear system it may be difficult to interpret the messages to get an idea of how long solve will take. If you have a polynomial system then there are other algorithms in Maple 10 that can be used. What is the input ?

For your example: LinearAlgebra:-LinearSolve(A.A, E); In general (solve Ax=B): LinearAlgebra:-LinearSolve(A,B); Keep in mind, this isn't going to work well beyond very small sizes, because the determinant of A (which appears in the denominator of the solution) blows up. If A is a 10x10 matrix with distinct entries, the determinant will be a polynomial with 10! =~ 3.6 million terms.

mpl files are typically just text files containing Maple source code. If the code has errors then you will get an error when you read it into Maple. You should not use Maple to edit mpl files. Use a text editor instead.

You might try some of the following:

f := 2*x + 3*y + z;
L := [op](f);   # list (indexed from 1)
R := rtable([op](f));  # rtable (array) indexed from 1
S := rtable(0..nops(f)-1, [op](f));  # rtable indexed from 0
T := table([op](f));  # hash table with T[i] = ith term

[op](f) might be slightly more efficient than [op(f)] for large f. I wouldn't use the lowercase array() because the underlying data structure is just a hash table with the overhead of bounds checking. rtables are usually more appropriate. If you need the coefficients and monomials in separate lists use the coeffs command.

To use the NAG routines, specify datatype=float[8]. The compiled routines are very fast, however I always seem to run out of precision, at which point Maple switches to software floats and becomes quite slow. Assuming you want an exact result, here are a few "tips": 1) storage=sparse is usually slower than the default rectangular storage. This is because no sparse algorithms were implemented so what you get is the artless strategy of a dense algorithm with the overhead of a sparse data structure. The data structure is actually a hash table, so I guess technically it supplies sparse storage, but you can't really do any operations that are not O(|A|). 2) The algorithms in LinearAlgebra do not pivot for sparsity or possibly even for size. It is recommended for performance that you convert your Matrix into a 'matrix' (ie: convert(A, 'matrix'); and use the old linalg package which was superceded by LinearAlgebra about 7 years ago. The command is "rref", for example: linalg[rref](A); 3) If you just want to solve a linear system, I suggest converting it to a set of equations ( 2*x[1] + 3*x[3] = 5, ...} and passing it to the solve command. This will call the fastest code yet, which is now over 20 years old. For matrices as small as yours, it will probably suffice. The sad truth is there has been little progress for doing exact sparse linear algebra in Maple for a while now. It's not that Maple can't do it, it's just that building something for the top level is a massive undertaking. I have some code which can probably solve your problems quite quickly. Let me know by private message if you need it.