Let's say that the matrices are named A and B. In the first worksheet, do

save A, B, "MyFile.m":

The MyFile can be anything, but the extension should probably be .m if the matrices are large, as that causes an encoding into a Maple-specific format. In the second worksheet, do

read "Myfile.m": 

At this point, the names A and B will be defined in the second worksheet to whatever they were in the first worksheet at the time of the save.

Because of the special evaluation rules of seq, add, and mul, you have to do

eval((f@Seq)(i, i= 1..2), Seq= seq);

and likewise for the add. The Seq here is just a dummy; it could be any name. And note that i does not get the usual protection against global interference that those special evaluation rules provide.

Here's a completely different method based on the relation being given by a boolean procedure of two arguments rather than by a set of ordered pairs. That is, for x, y in S, R(x,y) must return true if x is related to y and false otherwise.

Matrix methods are used to check symmetry and transitivity.

EqvRel:= (R::procedure, S::set)->
   ((A::Matrix)->
      (
         'reflexive' = andmap(x-> R(x,x), S),
         'symmetric' = ArrayTools:-IsZero(A - A^+),
         'transitive' = ArrayTools:-IsZero(A - signum~(A^2 + A))
      )
   )
   (Matrix(nops(S)$2, (i,j)-> `if`(R(S[i],S[j]), 1, 0), datatype= integer))
:

Examples:

Rel1:= (x,y)-> irem(x-y, 3) = 0:
S:= {'rand()' $ 30}:

EqvRel(Rel1, S);

     reflexive = true, symmetric = true, transitive = true

Rel2:= (x,y)-> x <= y:
EqvRel(Rel2, S);

     reflexive = true, symmetric = false, transitive = true

Rel3:= (x,y)-> x <> y:
EqvRel(Rel3, S);

    reflexive = false, symmetric = true, transitive = false

Your input array([1 1 1 2]) needs to be array([1,1,1,2]). Commas are always needed as separators on input.

The option maxmesh only applies to BVPs; this is an IVP.

The problem is that alias only creates rather-superficial pointers. They are intended mostly to de-clutter the display of expressions. They don't work programatically. In other words, in the following commands, diff never sees the alias A, it only sees f(x):

alias(A=f(x)):
diff(A^2, A);
Error, invalid input: diff received f(x), which is not valid for its 2nd argument
 

You can create more-substantial pointers with freeze, and dereference them with thaw:

restart:
A:= freeze(f(x)):
thaw(diff(A^2, A));
     2*f(x)

If you also want to use an alias for f(x) so that you never see it in output because you consider it private, that'd work fine.

You mentioned pointers to vectors. I don't see any vectors in your code. You also load LinearAlgebra, but you don't use anything from that package.

Here is a solution for the first problem. I solve the ODE as a parametrized IVP with parameter e and then fsolve for some values of e such that the right-side boundary condition is satisfied. You'll need to try some fiddling around with fsolve's options to get more values. Please let me know if there are any details below that you can't figure out.

Update: See a Reply below for vastly improved version of this worksheet where I show how to fsolve for all the eigenvalues in a given interval.

Download ShootingForEigenvalues.mw

Every closed hyperrectangle (henceforth: rectangle) in R^d with edges parallel to the coordinate axes can be uniquely specifed by what I'll call its principal diagonal, which is the unique diagonal with endpoints A and B such A[i] < B[i] for i= 1..d. A rectangle is to be divided into N^d congruent subrectangles. These can be indexed by i, 0 <= i < N^d. The following procedure returns the ith subrectangle:

alpha:= (N::posint, d::posint)-> 
   (i::nonnegint)-> convert(N^d+i, 'base', N)[..d]
:
Subrectangle:= proc(i::nonnegint, R::Rectangle(list,list), N::posint)
local A:= op(1,R), Alpha:= alpha(N,nops(A))(i), v, a;
   v:= (op(2,R) -~ A)/~N;
   a:= A +~ Alpha*~v;
   'Rectangle'(a, a +~ v)
end proc:

Subrectangle(5, Rectangle([0,0,0], [1,1,1]), 3);
     Rectangle([2/3, 1/3, 0], [1, 2/3, 1/3])

One possible such command is codegen[cost], but in this case the exponentiations (with small integer exponents) would be counted as an appropriate number of multiplications, because that's how the expression would be numerically evaluated.

There's no need to use arrow expressions; they just complicate things in this case. The following code does what you want. It works for any number of input expressions. It doesn't rely on any specific names being used such as x, y, or t.

linearise:= proc(S::list(algebraic), P::list(anyfunc(name)=algebraic))
local F:= lhs~(P), t:= op~({F[]})[];
   if nops([t]) > 1 then 
      error "Expected exactly 1 independent variable but received", t
   end if;
   diff~(F,t) =~ thaw(map(mtaylor, subs(F=~ freeze~(F), [S,P])[], 2))
end proc:

#Example usage
a:= x(t)*(y(t)-1);
b:= 4 - y(t)^2 - x(t)^2;
linearise([a,b], [x(t)= 0, y(t)= -2]);

Here's another procedure. This'll work in any number of dimensions, as long as the hyperfaces are parallel to the coordinate hyperplanes.

Vertices:= proc(P::list, Q::list) 
local x,p;
   seq([seq(x, x= p)], p= Iterator:-CartesianProduct(zip(`[]`,P,Q)[]))
end proc:

Vertices([1,1,2], [3,4,5]);

I don't know whether it can be done with geom3d.

I suggest that you stop using ancient and deprecated stats package and instead use the newer Statistics package. Then

x:= Statistics:-Sample(Uniform(0,1), 10);
add(x)^a;

Since the above is hardly worth doing, I think that what you actually want is

add(x_i^a, x_i= x);

Right-click on the Matrix or Vector output and select Browse.

The following does exactly what you're trying to do using a different ODE (because you didn't supply yours) and using a list of three operators as the first argument to plot3d. Because of the option remember, td is only called once (instead of three times) for each mesh point of the plot.

ans:= dsolve({diff(y(t),t$2)=y(t), y(0)=a, D(y)(0)=0}, numeric, parameters= [a]):
td:= proc(x,tt)
option remember;
   ans(parameters= [a= x]);
   eval([t, y(t), diff(y(t),t)], ans(tt))
end proc:

TD:= (k::{1,2,3})-> (x,tt)-> td(x,tt)[k]:
plot3d(TD~([1,2,3]), 0..3, 0..3, labels= [t, y, `y'`]);

In Maple, several binary infix operators {=, <, .., ^, ::, etc.} always have exactly two operands. For =, each of the operands can be anything---including NULL, a sequence, or another equation---as long as it's enclosed in parentheses. So (A = B) = C and A = (B = C) are both allowed. I realize that the display of the parentheses in printed output is not ideal.

This is implemented in the kernel, so it works regardless of which GUI (Standard, Classic, command-line) or GUI mode (document, extended typesetting, etc.) you use.

P:= proc(u::{function, indexed}, shift::anything:= 1)
local pos:= `if`(procname::indexed(posint), op(procname), 1);
   subsop(pos= op(pos, u) + shift, u)
end proc:

P[1](u(i,j));
                          u(i + 1, j)
P[2](u[i,j,k], 3);
                          u[i, j + 3, k]

This allows you to have a single operator that does it all. It doesn't care whether you've used square brackets or parentheses. It doesn't care how many indices there are. And you can change the shift amount from the default of 1,