restart:
  with(geometry):
  _EnvHorizontalName:=x:
  _EnvVerticalName:=y:
  x0 := (3+1/2)/2:
  x1 := x0: y1 := x0 + (1+1/4):
  point(P1, x1, y1):
  circle(C2, [P1, x0]):
  x3 := (3+1/2)/2-(3/4+20/1000)/2-1/2:
  y3 := 1+1/4: evalf([x3, y3]):
  point(P2, x3, 0):
  point(P3, x3, y3):
  line(L1, [P2, P3]):
  intersection(I1, L1, C2, [IP1, IP2]):
  tr := 1/8:
  x4 := x3 - tr:
  point(P4, x4, y4):
#
# Calculate all four possible solution circles
#
  circle(C4a, [point(P4a, x4, eval(y4, solve({distance(P4, P1)=tr+x1, y4<y1}))), tr]):
  circle(C4b, [point(P4b, x4, eval(y4, solve({distance(P4, P1)=x1-tr, y4<y1}))), tr]):
  circle(C4c, [point(P4c, x4, eval(y4, solve({distance(P4, P1)=tr+x1, y4>y1}))), tr]):
  circle(C4d, [point(P4d, x4, eval(y4, solve({distance(P4, P1)=x1-tr, y4>y1}))), tr]):
#
# Plot of all solutions
#
  draw([L1(color=blue), C2(color=green), C4a(color=red), C4b(color=red), C4c(color=red), C4d(color=red)]);
#
# Details of "lower" two solution circles
#
  draw([L1(color=blue), C2(color=green), C4a(color=red), C4b(color=red)], view=[0..1,1..2]);
  detail(C4a);
  detail(C4b);
#
# Details of "upper" two solution circles
#
  draw([L1(color=blue), C2(color=green), C4c(color=red), C4d(color=red)], view=[0..1,4..5]);
  detail(C4c);
  detail(C4d);

  restart;
#
# Define three signals as simple lists
#
  A:=[seq( [false$4, true$4][], j=1..4)]:
  B:=[seq( [false$2, true$2][], j=1..8)]:
  C:=[seq( [false$1, true$1][], j=1..16)]:
  myXOR:= (U,V)-> zip( `xor`, U ,V):
  toPW:=M->piecewise(seq( [t<j, `if`(M[j],1,0)][], j=1..numelems(M))):
#
# Plot relevant waveforms, offset in the y-direction
# by 2, and adjust the y-axis tickmarks
#
  plot( [ toPW(myXOR( myXOR(A, B), C) ),
          toPW(A)+2,
          toPW(B)+4,
          toPW(C)+6
        ],
        t=0..16,
        tickmarks=[ default,
                    [seq( j=j mod 2, j=0..7)]
                  ]
    );

  restart;
#
# Define three signals as simple lists
#
  A:=[seq( [0$4, 1$4][], j=1..4)]:
  B:=[seq( [0$2, 1$2][], j=1..8)]:
  C:=[seq( [0$1, 1$1][], j=1..16)]:
#
# Define a function which performs 'xor'
# on the elements of 2 lists
#
  myXOR:= (V,U) -> zip((i,j)-> (i+j) mod 2, V, U):
#
# Define a function which converts a
# list to a piecewise "waveform"
#
  toPW:=M->piecewise(seq( [t<j, M[j]][], j=1..numelems(M))):
#
# Plot relevant waveforms, offset in the y-direction
# by 2, and adjust the y-axis tickmarks
#
  plot( [ toPW(myXOR( myXOR(A, B), C) ),
          toPW(A)+2,
          toPW(B)+4,
          toPW(C)+6
        ],
        t=0..16,
        tickmarks=[ default,
                    [seq( j=j mod 2, j=0..7)]
                  ]
    );
 

restart;
with(VectorCalculus):
#
# Define three parameterized vectors in
# Cartesian coordinates (the default)
#
  p:=Vector(3, [ p__1(t), p__2(t), p__3(t) ]);
  q:=Vector(3, [ q__1(t), q__2(t), q__3(t) ]);
  w:=Vector(3, [ w__1(t), w__2(t), w__3(t) ]);
#
# Construct the integrand. NB this results in
# a scalar quantity
#
  IG:=(w.w +m^2)^(1/4)*(p.p+mu*2)^(1/4)/((w.w+a^2)^2*(p.p+a^2)*(q.q+b__2^2)^2*((p-q+w).(p-q+w)+b__1^2) );
#
# Set up the path integral along the vector 'p'
# OP seems to want to integrate along the vector
# from -Infinity to Infinity, so that the range
# of the parameterisation variable will have
# to ensure this - the following just uses
#
#            t=-Infinity..Infinity
#
# more or less as a placehoder
#
  PathInt(IG, [x,y,z]=Path(p, t=-Infinity..Infinity), 'inert');
#
# There appears to be no reason why these path
# integrals could not be nested
#

  restart;
  with(LinearAlgebra):
  with(plots):
  with(plottools):
#
# Set up matrix size and random number
# generators. Things start to slow down
# down a bit with more than 10*10 matrices,
# so be careful what you ask for!
#
  nR:=10:
  nC:=10:
  rv:=rand(0..1):
  s:=()->(-1)^rv():
  C:=RandomMatrix(nR, nC, generator=s):
  rR:=rand(1..nR):
  rC:=rand(1..nC):
#
# Define the procedure which changes a matrix
# entry, looking for a "lower energy"
# configuration
#
  doIter:=proc( mat::Matrix, rv1::procedure, rv2:=procedure)
                local augment, getE, oldE, newE,
                      matC:=mat, roll1, roll2,
                      count:=0:
              #
              # Define a procedure which "augments" the input
              # matrix so as to mimic periodic boundary conditions
              #
                augment:=proc( M::Matrix )
                               local augM;
                               augM:=Matrix( op([1,1],M)+2,
                                             op([1,2],M)+2,
                                             fill=0
                                           );
                               augM[2..op([1,1],M)+1, 2..op([1,2],M)+1]:=M;
                               augM[1, 2..op([1,2],M)+1]:= M[op([1,1],M),..]:
                               augM[op([1,1],M)+2, 2..op([1,2],M)+1]:= M[1,..]:
                               augM[2..op([1,1],M)+1,1]:=M[..,op([1,2],M)]:
                               augM[2..op([1,1],M)+1, op([1,2],M)+2]:=M[..,1]:
                               return augM
                        end proc;
              #
              # Define a simple energy function for the
              # augmented matrix
              #
                getE:= proc( M::Matrix)
                             add
                             ( add
                               ( (M[i-1,j]+M[i+1,j]+M[i,j-1]+M[i,j+1])*M[i,j],
                                 i=2..op([1,1],M)-1
                               ),
                               j=2..op([1,2],M)-1
                             )
                             +
                             4*(op([1,1],M)-2)*(op([1,2],M)-2);
                       end proc:
              #
              # Compute the "energy" of input matrix
              #
                oldE:=getE(augment(matC));
              #
              # Randomly change one entry in the input matrix
              # and check the new "energy". If it is <= to
              # the old "energy", then accept this new matrix
              # and exit.
              #
              # Restrict the number of attempts to 1000, just
              # to stop this getting stuck *permanently*
              #
                while true do
                      roll1:=rv1():
                      roll2:=rv2():
                      matC[roll1, roll2]:=-1*matC[roll1, roll2]:
                      newE:=getE(augment(matC)):
                      if   newE <= oldE
                      then return [newE, matC]
                      else matC[roll1, roll2]:=-1*matC[roll1, roll2]:
                      fi:
                      count:=count+1:
                      if count=1000
                      then return "ERROR"
                      fi:
              od:
       end proc:
#
# Define a plotting function, which I wouldn't have to
# do if the HeatMap() command could be used in animation,
# (and I still have no idea why it can't
#
  colors:=[red,blue]:
  gP:= M-> display
           ( [ seq
               ( seq
                 ( rectangle
                   ( [i-1,j-1],[i,j],
                     color=colors[(M[i,j]+3)/2]
                   ),
                   i=1..op([1,1],M)
                 ),
                 j=1..op([1,2],M)
               )
             ],
             scaling=constrained
           ):
#
# Initialise the loop and produce the starting
# plot
#
  p[1]:=C:
  q[1]:=gP(p[1]):
#
# Perform iterations. If doIter() is successful
# (ie produces a matrix) then generate a plot.
# If it is unsuccessful, then it will return a
# string, at which point this loop will terminate
#
# For a starting 10*10 Matrix, it takes ~150
# iterations through this loop to reach the
# ultimate "minimum energy" state, although
# it is not obvius to me that the final
# checkerboard (minimum energy) state would
# ever (necessarily) be obtained by the random
# process - getting stuck in some kind of
# cycle must (surely?) be a possibility
#
# So the upper limit on this loop is just to
# avoid excessive run times - you can raise
# it, but do so carefully
#
  for k from 2 by 1 to 200 do
      ans:= doIter(copy(p[k-1],deep), rR, rC);
      if   type( ans[2], Matrix )
      then p[k]:= ans[2]:
           E[k]:= ans[1]:
           q[k]:=gP(p[k]):
      else break
      fi
  od:
  printf("Actually did %d iterations\n",k-1):

Actually did 149 iterations

#
# Show the total "energy". This ought to be a
# decreasing function, but may or may not reach
# zero, depending on whether the lowext enrgy state
# has been achieved
#
  listplot( convert(E, list) );
#
# Animate the results which have been obtained.
# This can be a bit slow, depending on the size
# of the matrix, and the number of iteration
# steps which have been used. For 10*10 matrices
# and ~200 iterations, I would call it acceptable
# (ie ~30secs)
#
  display( convert(q, list), insequence=true);

restart;
with(ListTools):
M:=Matrix(4, 4, [ [ 1, 2, 3, 4 ],
                  [ 2, 3, 4, 5 ],
                  [ 3, 4, 5, 6 ],
                  [ 4, 5, 6, 7 ]
                ]
         ):
Ent:=[entries(M, nolist)]:
Ind:=[indices(M)]:
[{seq( [{j, seq( Ind[k], k in SearchAll(j, Ent))}[]], j in Ent)}[]];

  restart;
  interface(rtablesize=20):

#
# Generate a *list* each of whose entries is
# a 2x2 matrix
#
  A:= [ seq
        ( Matrix
          ( [ [ k,   k^2],
              [ 2*k, 2*k^2]
            ]
          ),
          k=1..10
        )
      ];
  whattype(A);

#
# Generate a *set* each of whose entries is
# a 2x2 matrix
#
  A:= { seq
        ( Matrix
          ( [ [ k,   k^2],
              [ 2*k, 2*k^2]
            ]
          ),
          k=1..10
        )
      };
  whattype(A);

#
# Generate a *table* each of whose entries is
# a 2x2 matrix
#
  A:= table
      ( [ seq
          ( Matrix
            ( [ [ k,   k^2],
                [ 2*k, 2*k^2]
              ]
            ),
            k=1..10
          )
        ]
      );
  whattype(op(A));

#
# Generate an *Array*, each of whose entries is
# a 2x2 Matrix
#
  A:= Array
      ( 1..10,
        k->Matrix( [ [ k,   k^2],
                     [ 2*k, 2*k^2]
                   ]
                 )
      );
  whattype(A);

#
# Generate a column Vector, each of whose entries
# is a 2x2 Matrix
#
  A:= Vector[column]
      ( 10,
        k->Matrix( [ [ k,   k^2],
                     [ 2*k, 2*k^2]
                   ]
                 )
      );
  whattype(A);

#
# Generate a row Vector, each of whose entries
# is a 2x2 Matrix
#
  A:= Vector[row]
      ( 10,
        k->Matrix( [ [ k,   k^2],
                     [ 2*k, 2*k^2]
                   ]
                 )
      );
  whattype(A);

#
# Generate a matrix, whose diagonal entries
# are the required 2x2 Matrices
#
  A:= Matrix( 10,
              10,
              (i,j)->`if`( i=j,
                           Matrix( [ [ i,   i^2],
                                    [ 2*i, 2*i^2]
                                   ]
                                 ),
                           0
                         )
            );
  whattype(A);