What does x | i mean? It's part of the argument of JacobiTheta0.

@testht06 The identity matrix in Maple has no predefined symbolic representation; it needs to be constructed to whatever size is appropriate in context. Capital I is the imaginary unit rather than a matrix.

This works:

A:= M + LinearAlgebra:-IdentityMatrix(op([1,1], M)):
Gausselim(A, 'rank') mod 2:
rank;

@tomleslie A restart command needs to be in its own execution group to be fully effective. When it's in an execution group with other commands, the results are unpredictable. This is not clearly documented at ?restart, which merely says that the restart must be on its own "line". If you put the restart in its own execution group, my counterexample will be manifested.

@testht06 Yes, I know how to perform division in finite fields. That doesn't change the fact that the sum of a row could be zero. 

@sipa Please repost this as a separate Question, a new thread. The syntax for plotting ODE solutions is very different from that for plotting PDE solutions. There's no way to adapt the code that I give above to an ODE situation.

For your item (3), do you want the division to be in the finite field? I'm not sure if that's always possible. Couldn't the sum squared be 0 in the field?

I notice that you index your matrix columns starting with 0 rather than 1. But matrices (as opposed to arrays) are always indexed from 1. Is that going to cause a problem for you?

Since, as Tom points out, this can't really be done, I'd like to ask what your ultimate goal is. There's probably some other way to achieve it.

@tomleslie 

Your solution relies on the terms of the polynomial being put in the same order by both op and coeffs, but it doesn't necessarily happen that way. Simple example:

restart:
p:= y+2*x;
op(p);
coeffs(p);

@acer Why did you use assign? Why not just

UpdatingRead:= proc(fn::string) ... end proc?

@mela 

The following procedure, MostConstrainedVariable, will get you started towards automating this process. If fsolve actually returns a solution, and you specified ranges for all the variables, this procedure will tell which variable is closest (in relative terms) to violating its range. The returned value is the variable, its value, its range, and the fraction of that range from the value to its nearest boundary. Remember, this procedure only works for cases where fsolve does return a solution.

Download MostConstrainedVariable.mw

Also, I'd like to point out that Maple's fsolve doesn't require you to specify ranges for the variables. You can specify initial guesses just like you say that you do in Matlab. Or you can specify both a range and an initial guess for a variable. Or you can specify nothing, although I don't advise doing that for a system as complicated as yours.

@aamirkhan Yes, simplify is designed to be exact. In other words, applying simplify to an expression results in another expression which is (or is supposed to be) mathematically identical to the first (numerically equal over all possible variable assignments over the complex numbers). If your expression contains floating-point numbers, that equality may be modulo some round-off error. The chance that that round-off error is manifested as a dramatic reduction in the size of the expression is infinitesimal.

@aamirkhan 

For some strange reason, fnormal is designed to automatically map over lists, sets, etc., but not over Matrices, Vectors, etc. Your M1f is a Matrix. We can correct for that deficiency in fnormal by using map like this:

M1ff:= map(fnormal, M1f, Digits, 5e-7);

where the 5e-7 indicates that any value with magnitude less than 5e-7 is to be changed to 0. 

You can correct fnormal with overload so that it automatically maps over Matrices, etc., like this

unprotect(fnormal):
fnormal:= overload([
     proc(e::rtable) option overload; map(fnormal, e, _rest) end proc,
     fnormal:-ModuleApply
]):
protect(fnormal):

Then proceed to use fnormal as it is documented at ?fnormal except that now the first argument can also be a Matrix, Vector, etc.:

M1ff:= fnormal(M1f, Digits, 5e-7);

You say that you want t[2] on the z-axis. I see no t[2] in your function.

You need to correct T[g/T[0]] to T[g]/T[0].

Please post a worksheet showing the command that you type and the resulting plot. One example should be sufficient, but include as many as you wish. Use the green uparrow (last item on the second row of the toolbar in the MaplePrimes editor) to upload and display the worksheet directly in your Reply.

@MDD I think so, but I am not sure. You can run debug on any procedure or group of procedures from Groebner. Then run an eliminate command on polynomials. If any of the Groebner procedures that you listed are called, you'll be notified of thatby the debug.