I got it to work using Groebner:-Basis instead of Primfield. I assume that you still want the fourth root of the starting matrix; that's what's computed in the worksheet below. Other roots could be computed just as well.

I couldn't have done this without the crucial idea supplied by Christian Wolinski, so I promoted his Reply with that idea to an Answer so that I could give it a vote up.

Sorry, for some reason this worksheet won't display inline, but here's the link:

Matrix_powers_finite_field_2.mw

It seems not possible to get a solution for some negative values of T. The following code will produce a matrix of the solutions that it is possible to get.

restart:
Digits:= 15:

eq1 := PI_T = alpha*M*(kappa/(1+kappa-sigma)-1):
eq2 := l = alpha*((sigma-1)*kappa/(1+kappa-sigma)+1):
eq3 := M = phi_c^(-kappa)*F:
eq4 := e*F = PI_T*v+T:
eq5 := M*l+L_A = L_s:
eq6 := F*e+A = A_s:
eq7 := A = (1-beta)*(L_s+(1-v)*PI_T-T):
eq8 := A_s = L_A:
eq9 := p = sigma/((sigma-1)*phi):
eq10 := P_Y = M^(lambda/(1-sigma))*p:
eq11 := L_s = N*(theta*P_Y^beta)^(-1/delta):
eq12 := phi = (kappa/(1+kappa-sigma))^(1/(sigma-1))*phi_c:
eq13 := U = delta*theta*L_s/((1+delta)*P_Y^beta)+((1-v)*PI_T-T)/P_Y^beta:

Params:= {N = 1, v = 1, beta = .75, lambda = 1, sigma = 5, kappa = 4.8, delta = 2, theta = 2.591350635, alpha = 0.3998699153e-1, e = 0.8333333332e-1}:

Init_Values:= {A_s = .2500000000, l = .9996747882, PI_T = 0.8333333332e-1, L_A = .2500000000, phi = 1.878101496, M = .4168022157, A = .1666666667, p = .6655657336, P_Y = .8283396488, L_s = 2/3, phi_c = 1.2, F = 1, U = 1.326439132 }:

V:= <T, A_s, l, PI_T, L_A, phi, M, A, p, P_Y, L_s, phi_c, F, U>^+:

SOLs:= table():
for _T from -0.3 by 0.01 to 0.3 do     
     s:= fsolve(eval({eq||(1..13)}, Params union {T = _T}), Init_Values);
     if op(0, eval(s, 1)) = fsolve then next end if;
     SOLs[_T]:= eval(V, s union {T= _T})
end do:

sol_matrix:= Matrix(<entries(SOLs, 'nolist', 'indexorder')>, datatype= float[8]);

simplify(P(epsilon), {epsilon^2= 0});

This will change epsilon^2 and all higher powers of epsilon to 0.

Rather than a[k]:= 'a[k]', use unassign('a[k]').

If you're going to remove a lot of entries from a table, it may be better to create a new table instead:

a:= table(remove(e-> lhs(e)::even, op(2, eval(a))));

plots:-shadebetween(4, x^2, x= -2..2);

A subs command requires at least three parts: A, B, C:

subs(A= B, C);

You're missing the A, which can be any of K[4], W(x), omega, etc. I suspect that you intend for A to be K[4], but I'm not sure. If this is what you intend, then

EquConFour:= subs(
     K[4]= int(W(kappa)*exp^(-J*kappa*x), kappa= -infinity..infinity)/2/Pi,
     EquCon
)

If you want a single selection, then do

rand(0..100)();

If you want multiple selections, then do

R:= rand(0..100):

and call R() every time you want a number.

The algorithms for implicitplot and implicitplot3d are different. The latter is purely numeric and can't consider things like factorization. The only accuracy control that you have is the grid or numpoints option.

Define this procedure:

RemoveQuotes:= (e::uneval)-> subsindets(e, uneval, eval, 1):

Then copy-and-paste the output of the tutor (excluding the terminal semicolon) as the input of that procedure:

RemoveQuotes( paste here );

Then copy-and-paste the resulting output as the input at the next prompt. Then edit away to your heart's content.

I find define very difficult to use and prefer to write a procedure directly. Try this:

E:= proc(e::algebraic)
local a,b;
     if not hastype(e, RV) then e
     elif e::RV then 'procname'(e)
     elif e::`+` then map(thisproc, e)
     elif e::`*` then
          (a,b):= selectremove(hastype, e, RV);
          b*thisproc(expand(a))
     elif e::`^` and op(2,e)::posint then thisproc(expand(e))
     else 'procname'(e)
     end if
end proc;

[Code of E edited since original posting. See Replies.]

By the way, your type RV can be corrected and simplified to

TypeTools:-AddType(
     RV,
     {RandomVariable,
       'RandomVariable^posint',
        '`*`'({RandomVariable, 'RandomVariable^posint'})
     }
);

with no need to use the module prefix Statistics:-, regardless of whether the package has been loaded. The correction is that it now allow products of powers of random variables, and not every expression of type `*` has two operands.

According to the standard definitions (see Wikipedia) a walk is not a path (but a path is a walk). However, I think that an animation highlighting both the vertices and the edges being travelled is easier to understand.

restart:

#Kitonum's procedure, slightly modified.
RandomWalkOnGraph:= proc(G::Graph, s::{symbol,integer,string}, N::integer)
local S, i;
uses GraphTheory, RandomTools;
     S[1]:= s;
     for i from 2 to N do  S[i]:= Generate(choose(Neighbors(G, S[i-1])))  end do;
     convert(S, list)
end proc:

G:= GraphTheory[Graph](GraphTheory[Trail](1,2,3,4,5,6,4,7,8,2)):
W:= RandomWalkOnGraph(G,2,20);

     W := [2, 1, 2, 1, 2, 8, 7, 4, 5, 4, 5, 4, 7, 4, 6, 5, 4, 6, 4, 6]

AnimateWalk:= proc(G::Graph, W::list({symbol,integer,string}))
uses GT= GraphTheory;
local k, n:= nops(W)-1, P:= Vector(3*n), G1;
     for k to n do
          G1:= GT:-Graph(GT:-Edges(G));
          GT:-HighlightVertex(G1, W[k], red);
          P[3*k-2]:= GT:-DrawGraph(G1);
          GT:-HighlightEdges(G1, {W[k],W[k+1]});
          P[3*k-1]:= GT:-DrawGraph(G1);
          GT:-HighlightVertex(G1, W[k+1], red);
          P[3*k]:= GT:-DrawGraph(G1)
     end do;
     plots:-display(convert(P,list), insequence)
end proc:

AnimateWalk(G,W);

This will apply the rule to all powered parameters. (Sorry, this isn't elegant.)

thaw(
     subsindets[2](subsindets(z2, {identical(alpha), identical(alpha)^anything}, freeze), `^`, 1, op)
);

If you would like to apply the rule only to parameters named a[something], let me know.

Assuming that the sizes of the component matrices are compatible, the following simple procedure will flatten (or resolve) it into an ordinary Matrix.

Flatten:= proc(M::Matrix(Matrix))
local i,j;
     `<|>`(seq(<seq(M[i,j], i= 1..op([1,1], M))>, j= 1..op([1,2], M)))
end proc:

B, Flatten(B);

In addition to Acer's great suggestions, there's another tool that sometimes can help you track down where the memory is being used: kernelopts(memusage). This breaks the memory usage down into 63 categories and gives the number of items in each category and the total memory usage for the category. You can use subtraction to compare this matrix before and after a single call to MyProc.

Regarding Profiling: I'm a big fan of it, but I think that profiling a single suspect procedure is worthless if that procedure calls any other procedure, including Maple library procedures. Profile everything with CodeTools:-Profiling:-Profile();

In the differential equations, change every occurence of x to x(t) and every occurence of y to y(t). You should get this plot: