You show "Return", which is neither return nor RETURN. There is no Return.

You wrote that you wanted it nonsingular. So why are you returning it if the determinant is zero? Isn't that backwards?

Also, what is it returning from? Where's the proc? Did you want something recursive, perhaps like this?

 (**) MakeKeyMatrix:=proc(n)
 (**) local KeyMatrix;
 (**) uses LinearAlgebra:-Modular,LinearAlgebra;
 (**) KeyMatrix:=Mod(127,RandomMatrix(n),integer);
 (**) if Determinant(127,KeyMatrix)=0 then
 (**)   return MakeKeyMatrix(n);
 (**) else
 (**)   return KeyMatrix;
 (**) end if;
 (**) end proc:

It seems likely that what you are (repeatedly) assigning to t are actually sequences, not lists. A list can be constructed like a sequence, but with [ ] brackets.

It's not a great way to code, by repeatedly appending to the end of a sequence or list. But if you insist on doing it that way, then make sure that t really is a list (not a sequence) so that you can index into it when it has just one element.

I am guessing that you are trying something like this,

> t:=NULL;
                                     t :=

> t:=t,7;
                                    t := 7

> t[1]; # oops
                                     7[1]

> t:=t,13;
                                  t := 7, 13

> t[2];
                                      13

This is not efficient or pretty, but may be a close modification of your technique.

> restart:

> t:=NULL;
                                     t :=

> t:=t,7;
                                    t := 7

> [t][1];
                                       7

> t:=t,13;
                                   t := 7, 13

> [t][2];
                                       13

Alternatively,

> restart:
> t:=[NULL];
                                    t := []

> t := [op(t),7];
                                   t := [7]

> t[1];
                                       7

> t := [op(t),13];
                                  t := [7, 13]

> t[2];
                                       13

> Student[Calculus1][RiemannSum](f(x), x=a..b,method=midpoint,output=sum);
                              /  9                           \
                              |-----                         |
                              | \          (i + 1/2) (b - a) |
                 1/10 (b - a) |  )   f(a + -----------------)|
                              | /                 10         |
                              |-----                         |
                              \i = 0                         /

> S:=Student[Calculus1][RiemannSum](f(x), x=a..b,
>       method=midpoint,output=sum,partition=N);
                              /N - 1                         \
                              |-----                         |
                              | \          (i + 1/2) (b - a) |
                      (b - a) |  )   f(a + -----------------)|
                              | /                  N         |
                              |-----                         |
                              \i = 0                         /
                 S := ----------------------------------------
                                         N

> eval(S,[f=sin,a=0,b=1]);
                              N - 1
                              -----
                               \        i + 1/2
                                )   sin(-------)
                               /           N
                              -----
                              i = 0
                              ------------------
                                      N

> evalf(Limit(%,N=infinity));
                                 0.4596976941

> int(sin(x),x=0..1.0);
                                 0.4596976941

You can tell `int` not to take time looking for discontinuities.

> int(s,y=evalf(-Pi)..2.338412729,continuous);
                             4.931498685

Alternatively, you can specify use of the fundamental theorem of calculus as the method.

> int(s,y=evalf(-Pi)..2.338412729,method=FTOC);
                             4.931498685

Both of those are fast here, taking less than 1 second on my machine.

I believe that the reason your initial attempt took so very long is that it was dispatching to evalf/Int, and that numerical routine was getting all caught up looking for discontinuities. (For definite integration with floating-point end-points, `int` passes off to evalf/Int.) That behaviour of `evalf/Int` to search for discontinuities of an expression can also be forcibly avoided by passing it a procedure instead of an expression.

> st:=time(): evalf(Int(Y->eval(s,y=Y),-Pi..2.338412729));
> time()-st;
                              4.931498685
                                 4.016

The code above shows that numerical quadrature is slower here even when not bogged down searching for discontinuities. Numerical quadrature can also be bypassed without having to force the 'method' or 'continuous' options of `int`. It can be done simply by passing exact end-points.

> st:=time(): Digits:=15:
>             int(s,y=-Pi..convert(2.338412729,rational,exact));
>             Digits:=10:
> time()-st;
                            4.9314986849254
                                 0.802

Note the capital M for Matrix.

c^%T;

LinearAlgebra:-Transpose(c);

linalg:-transpose(c);

Does it hang if you have something in your .mapleinit, or if you invoke it with some innocent optional command? Eg,

   maple -c "1:"

Is it possible that, when you entered the definition of hphi2 in 2D Math, you inadvertantly put a space between the cos and the left-bracket ( for the term cos(M*phi)?

> Az:=.15:
> if Az < .25 then print("O.K") else print("N.G") end if;
                                     "O.K"

> Az:=.15*Unit(km):
> if convert(convert(Az,units,mm),unit_free) < .25 then
>   print("O.K");
> else
>   print("N.G");
> end if;
                                     "N.G"

You could read the Intat help page, or hit it with value().

There are doubtless more efficient ways.

> M:=Matrix([[2, 5], [4, 0], [3, 9]]);
                                   [2    5]
                                   [      ]
                              M := [4    0]
                                   [      ]
                                   [3    9]

> F:=A->Matrix(op(1,A),
>    map([rhs],sort([seq(A[i,1]=A[i,1..-1],i=1..op([1,1],A))],
>                   (a,b)->lhs(a) < lhs(b)))):

> F(M);
                                   [2    5]
                                   [      ]
                                   [3    9]
                                   [      ]
                                   [4    0]

> Matrix([Vector(combinat:-randperm([seq(i,i=1..3)])),
>         LinearAlgebra:-RandomVector(3,generator=1..3)]);

                               [2    1]
                               [      ]
                               [3    2]
                               [      ]
                               [1    1]

Note that a name with assumptions on it is not of type global.

> restart:
> type(M,`global`);
                                     true

> assume(M::real):
> type(M,`global`);
                                     false

So, yes, there could be confusion for you there. If you upload your worksheet to this site (green up-arrow in the editing toolbar) then someone might be able to sort it for you.

Whenever you place new assumptions on a name Maple produces a new "version" of that name (like an escaped local). Instances of the name produced under different assumptions are not recognized as being the same.

> assume(M::real):
> x := 2*M;
                                   x := 2 M~

> additionally(M::positive);
> y := 6*M;
                                   y := 6 M~

> y+x;
                                  6 M~ + 2 M~

> has(x,M);
                                     false


> restart:
> assume(M::real):
> additionally(M::positive);
> x := 2*M;
                                   x := 2 M~

> y := 6*M;
                                   y := 6 M~

> y+x;
                                     8 M~

> has(x,M);
                                     true

Now, I'm not saying that you did things in that exact order. But it sounds as if at least some instances of `M` were assigned under different assumptions. You could try to find those and fix by ensuring that all such assignments are done with the very same set of assumptions. (Executing worksheet lines "out of order" is a typical way to get into such a minor muddle.)

The left-hand-side parameters of an operator (procedure parameters) should be of type symbol and not merely of type name (ie, index names).

By default, a 2D Math subscripted name is an indexed name. Ie, the subscripted name that you get by typing B_1 is actually B[1] underneath.

Here are two alternative ways around this.

1) Make the parameter names simply B and C, and do nothing to the B_1 and C_1 inside the operator.

2) Toggle all instances of B_1 in both the parameter list as well as the operator body to be "atomic identifiers" using the right-click context sensitive menu choice 2-D Math -> Convert To -> Atmomic Identifier. Select all of B_1 before right-clicking. Note that it would be just as easy to instead use two other non-subscripted formal parameter names instead, and get the same effect here.

Which solution you choose may depend on what you intend to pass to the operator. If you intend on passing it argument objects which are indexed strucures (Vectors, lists, tables) then go with 1). If you intend on passing scalar objects then go with  2).

A few minor comments on efficiency:

Rather than pass the very same add() invocation in two different arguments to NLPSolve, you might try assigning the result of the add() to a variable beforehand, and then passing that assigned variable. (add does not have option remember, AFAIK, so it is unnecessary duplication to have the same addition get computed more than once.)

This next may only matter if your numeric Matrices are very large. Take the Phi.C.f computation, and look at the number of arithmetic operations done. If Phi and C are Matrices and f is a Vector then if can make a measurable difference in performance if you compute Phi.(C.f) instead of (Phi.C).f.