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:

Here's my improved version of Kitonum's ContoursWithLabels. Now there's no specific variable name dependence, and I've added separate keyword options to control each plot.

Download ContoursWithLabels.mw

plots:-polygonplot3d([A1, A2, A3]);

Here's an easier way for your whole code:

block:= ListTools:-Reverse~(Bits:-Split~(convert(message, bytes), bits= 8)):
Block:= Matrix(block):
Block[..,3..7]:= 1 -~ Block[..,3..7]:
cblock:= convert(Block, listlist);

Depending on what your next operation is, it may be better to skip the last command and just keep it in Matrix form.

Here's a variant that uses the special ways that Matrices can be indexed:

SelectRealRows:= (M::Matrix)-> M[remove(i-> has(M[i,..], I), [$1..op([1,1],M)]), ..];

The reason is that filled is an option to the plot command, not an option to the PLOT data structure. In other words, plot does some computation to construct the polygon to fill; it's not the plot renderer that does it.

Here's a workaround. This should work for any 2D plot with a single curve. (But, mind you, this is quick-and-dirty; I only tested it on DensityPlots.)

FillPlot:= proc(P::specfunc(PLOT))
local
     C:= indets(P, specfunc(CURVES))[1],
     A:= convert(op(1,C), listlist)
;
     PLOT(
          POLYGONS(
               [[A[1][1],0], A[], [A[-1][1],0], [A[1,1],0]],
               LEGEND("__never_display_this_legend_entry"),
               STYLE(PATCHNOGRID),
               TRANSPARENCY(0.4),
               indets(C, specfunc(COLOUR))[]
          ),
          op(P)
     )
end proc:    

FillPlot(DensityPlot(X));       

To retain the symbolic constants, use solve instead of fsolve.

There is a way to create true symbolic constants, but you can't get fsolve to express solutions in terms of them; fsolve is strictly numeric.

Two ways. In the first, change g to

g:= proc(t)
local j;
     if not t::numeric then return 'procname'(t) end if;
    ... the rest of g ...

And change the plot command to

plot([g(t)[1], g(t)[2], t= 0..1]);

The second way is to simplify everything:

f0:= t-> (t, 3.9*t*(1-t)):
IFS:= (i,x,y)-> piecewise(i=0, [y,x], i=1, [x,y+1], i=2, [x,y]+~1, i=3, [2-y,1-x])[]/2:
g:= t-> IFS(trunc(4*t), f0(frac(4*t))):
plot([g(t)[1], g(t)[2], t= 0..1]);