Apply the command fnormal to the results of floating-point computation before making an equality comparison with zero. 

plots:-display([
     plot(p),
     plots:-textplot(
          zip((a,b)-> [a[], ""||b], p, [$1..nops(p)]),
          align= {above, left}
     )
]);

If you cancel theta[m]s in your x, you get

From this, it is clear that x > 1/4 for the range of phi that you are using. So that plot showing x going to 0 is bogus. (Note: x is the vertical axis---labelled "Cooperation level.")

Where are you reading the "description of the function" that you mention? Is it at ?IsSimilar? What version of Maple? Looking in Maple 16 help, I see no mention of a third argument P. My help says that the third argument can be output= 'C' (the C being literal). Thus, to get P, you can do

P:= LinearAlgebra:-IsSimilar(A, B, output= 'C');

No initialization of P is required.

nops does not work on tables. Use numelems instead.

This procedure is not exactly elegant, but it seems a difficult problem, and this procedure solves it in full generality. It will handle any expression with products that have any number of factors that are Sums. It generalizes from two to any number by applying subsindets recursively. It generates a unique global index name for every distinct Sum. The sequence of indices defaults to _j1, _j2, ..., and the default prefix _j change be changed by using the indexprefix option to the procedure.

ProductOfSumsToSumOfProducts:= proc(expr, {indexprefix::name:= :-_j})
global _indexsuffix;
local
     r, #return value

     #extracts the factors of a product that are Sums
     Extractor:= P-> select(type, {op(P)}, specfunc(anything, Sum)),

     #selects products having multiple Sum factors
     Selector:= And(`*`, satisfies(P-> nops(Extractor(P)) > 1)),

     #transforms a product containing two or more Sums by combining
     #an arbitrary two of them
     Transformer:= proc(P::`*`)
     local
          S:= Extractor(P), S1, S2, S12,
          j1:= indexprefix||(:-_indexsuffix+1), #1st index
          j2:= indexprefix||(:-_indexsuffix+2)  #2nd index
     ;
          :-_indexsuffix:= :-_indexsuffix+2;
          S1,S2:= S[1..2][]; #arbitrary 1st two
          #combine them
          S12:= Sum(
               Sum(
                    subs(op([2,1],S1)= j1, op(1,S1)) * subs(op([2,1],S2)= j2, op(1,S2)),
                    j1= op([2,2],S1)
               ), j2= op([2,2],S2)
          );
          algsubs(S1*S2= S12, P)
     end proc
;
     if assigned(:-_indexsuffix) then  unprotect(':-_indexsuffix')  else  :-_indexsuffix:= 0  end if;
     r:= expr;
     while indets(r, Selector) <> {} do  r:= subsindets(r, Selector, Transformer)  end do;
     protect(':-_indexsuffix');
     r
end proc:

Example of use:

e:= Sum(x[k],k=1..n)*Sum(y[k],k=1..m)*Sum(z[k],k=1..p);

ProductOfSumsToSumOfProducts(e, indexprefix= j);

Here's a standard technique that works in almost any language: Put your program in a while loop. After your program, but still in the loop, get the user's input regarding continuing. Like this:

UserProgram:= proc()
local Again:= Y;
     while Again in {Y,y} do

          # Your original program here
          # ...
          # end of your original program
       
          Again:= "";
          while not Again in {Y,y,N,n} do
               Again:= readstat(cat(Again, "Run again? (Y/N) "));
               if not Again in {Y,y,N,n} then  Again:= "Invalid input. " end if
          end do
     end do
end proc;

That being said, I find the dialog boxes produced by readstat (and the whole family of a commands for taking user input) very awkward. If you want to write Maple programs that process user input from the keyboard, use a Maplet. But that's like learning a whole new language.

If your user wants to run the program again, they can move the mouse pointer up to the program call, change the parameters, and press Enter.

Here are two ways that you can use plots:-display for this. The first prints the plot one after the other. The second uses the fancy Array form of display, and a display inside the display, to display the plots horizontally. The only library command that can do this is display.

Download Array_plot.mw

The first argument to animate needs to be the name of a plot command. Otherwise, how would it know that you want to animate a tube plot rather than, say, a space curve? The second argument is a list containing the arguments to that plot command, including the animation parameter. The third argument is a declaration of the animation parameter and its range. Then there are optional options, and options that apply to the plot as a whole.

plots:-animate(
     plots:-tubeplot,
     [[cos(t+phi), sin(t+phi), t], t= 0..6*Pi, numpoints=200],
     phi= 0..2*Pi,
     frames= 20,
     orientation= [40,70], shading= z, axes= box, style= patchnogrid
);

There are a vast array of gridline options described on the help pages ?plot,axis and ?plot,tickmarks. You can control the number of gridlines and their spacing (also color and linestyle). The spacing does not need to be even; you can place them at precise axis intersections. You can control both the major and minor gridlines.

Dumb example:

plot(0, -2..2, axis[1]= [gridlines= 20, color= red], axis[2]= [gridlines= 15, color= green]);

If after reading those help pages you still need help, please ask further.

There is no doubt a vast literature on this very important problem of finding the connected components of a graph. However, without researching any of that literature, I just now came up with this queue-based algorithm. So, there is good chance that there's a faster algorithm. But the below works. (There are also the existing commands GraphTheory:-ConnectedComponents and networks[components] that do the job, but I assumed that the goal was to write your own.)

You may well ask Since the order doesn't matter, why use a queue when you could just add/subtract items one at a time from a set? Good question. The answer is that the latter are very expensive operations in Maple; each add/subtract requires the reconstruction of the entire set. Queues (which are implemented as tables) work much faster.

I urge you not to use networks; it has been deprecated in favor of GraphTheory.

My procedure processes your table-based graph representation and returns a partition of its vertices into the connected components. I leave it to you to convert the blocks of this partition back into table-based graphs, which is trivial. If you have any trouble with that, let me know. My procedure can be easily modified to work with the GraphTheory representation of graphs. Let me know if you want that. Despite the title of your Question, my procedure is not limited to two components.

#Given a graph, return a partition of its vertices into the
#connected components.

ConnectedComponents:= proc(G::table)
local
     Q, #SimpleQueue of vertices to be visited in current component
     Comps:= table(), #set of components to return
     C,  #component currently being built
     v, #current vertex
     V:= {indices(G, 'nolist')}, #vertices remaining to be visited
     k #component index
;
     for k while V <> {} do
          #new component
          C:= table([V[1]=V[1]]);
          Q:= SimpleQueue(G[V[1]][]); #Queue the neighbors of new vertex.
          while not Q:-empty() do
               v:= Q:-dequeue();
               if assigned(C[v]) then next end if; #We already visited it.
               C[v]:= v; #Add this one to the current component.
               Q:-enqueue~(G[v]) #Queue all its neighbors.
          end do;
          #Convert to a set.
          C:= {indices(C, 'nolist')};
          V:= V minus C;
          Comps[k]:= C
     end do;
     {entries(Comps, 'nolist')}
end proc:      

ConnectedComponents(T);

Change the line

from i = 1 to N do

to

for i from 1 to N do

A dsolve(..., numeric) solution can be obtained very quickly. Just put the option numeric at the end of your dsolve command. Remove the evalf. To plot the solutions, do

plots:-odeplot(coz1, [x,X11(x)]);

and likewise for X22 and X33.

Let's define the error, e, as the absolute value of the difference between the true value and the computed approximation. Then for an iterative method we have a sequence of errors e[1], e[2], .... We say that the order of convergence is q if limit(e[n+1]/e[n]^q, n= infinity) approaches a positive constant u. (See the Wikipedia article "Rate of convergence".) It follows that for sufficiently large n, ln(e[n+1]) ~ ln(u) + q*ln(e[n]), where ~ denotes approximate equality. This is a linear equation in q, so we can use the e as the data for linear regression to approximate the "slope" q.

The warning is due to several of iterates having order of convergence exactly equal to 3 (as far as floating-point computation is concerned). The warning can be ignored.

Download Halley.mw

Welcome. I moved your Question from the Posts category to the Questions category.

Try typing this. Type it explicitly; do not use palettes or menus.

int(exp(2*x), x);