Any three non-collinear points uniquely determine a circle (and, of course, a triangle). (That existence and uniqueness of the circle was surprising to me when I first learned of it, so many decades ago.) So, we just need the center and radius of the circle determined by the vertices of the triangle. This can be done easily with the geometry package. Below, I generate three points at random, find the circle, and plot the circle, triangle, and circumcenter.

restart:
R:= rand(-9..9):
G:= geometry:
G:-circle(C1, G:-point~([A,B,C], '['R()'$3]'$2)):
Triang:= G:-coordinates~([A, B, C, A]);
         Triang := [[5, -8], [3, 2], [-4, -6], [5, -8]]

Ccen:= G:-coordinates(G:-center(C1));
                              [109  -305]
                      Ccen := [---, ----]
                              [86    86 ]
R:= G:-radius(C1);
                      1         (1/2)     (1/2)
                R := ---- 124865      3698     
                     3698                      
plot(
    [[(R*~[cos,sin](t) +~ Ccen)[], t= -Pi..Pi], Triang, [Ccen]], 
    style= [line$2, point], symbol= solidcircle, scaling= constrained,
    thickness= 5
);

If you want to work through (with Maple) the algebra that underlies the computations done by geometry above, let me know; it's not terribly complicated.

One way is like this:

x:= <1, 2, 3, 4>;
k:= Matrix(4, 4, (i,j)-> x[i]+x[j]);

Nested for loops could also be used, but it's not ideal.

Do you want the matrix to contain only the coefficients of the power-3 terms of the equations? In that case, do

KKNL:= Matrix((nops(EqML), nops(Var[4])), (i,j)-> coeff(EqML[i], Var[4][j]))

Acer's Answer is good. I just want to add what I think may have contributed to your original confusion about type indexed versus type symbol: When option symbol= a is used with a Matrix command, the nonzero entries of the resulting matrix are type indexed, not type symbol. It is the a (alone) which is the symbol; not the entries such as a[1,1].

Your desired form of the solution contains several factors that don't depend on t. These factors could be absorbed into the constants of integration without changing the correctness of the solution.

I agree with @tomleslie that you should use add instead of sum and local indices and shouldn't use unevaluation quotes. But I think the more-fundamental problem is using catenation rather than true indexing (such as A[i,j,k]). So, I'd do this (in 1D input):

p:= (i,j)-> local k; add(A[i,j,k], k= 1..3):
q:= i-> local j; add(p(i,j), j= 1..3):
r:= add(q(i), i= 1..3);

I can import and parse your 2800x2800 matrix of rational functions that you posted as file "f.txt" in only 7 seconds, like this:

restart:
M:= CodeTools:-Usage(
    subsindets(
        ImportMatrix("C:/Users/carlj/Downloads/f.txt"),
        string,
        parse
    )
):
memory used=185.79MiB, alloc change=195.61MiB, 
cpu time=6.69s, real time=6.18s, gc time=1.22s

upperbound(M);
                           2800, 2800

#Example entry:
M[9,9];
          /                                         
    1     |      2                                  
--------- |a[1] h  (K__ref Kuy0[1] + K__ref Kuyb[1])
 3        \                                         
h  K__ref                                           

                      3                                            
     144 a[1] K__ref h  #mover(mi("A"),mo("&uminus0;"))[6, 6, 1]   
   + ----------------------------------------------------------- - 
                                b[1]                               

        1       /     
  ------------- \a[1] 
  17 b[1] ro[1]       

                                                            2 
  #mover(mi("ω",fontstyle = "normal"),mo("&uminus0;"))  

          2        \
  K__ref h  I__0[1]/

                                                 \
     1        2                                  |
   + -- b[1] h  (K__ref Kux0[1] + K__ref Kuxa[1])|
     17                                          /

The MathML operators shown above (mover, mo, mi, etc.) won't appear in ordinary prettyprinted output. Their appearence above is an artifact of posting on MaplePrimes.

The subsindets(..., string, parse) is needed because symbolic expressions are being imported as strings, which isn't really a problem because parse very quickly converts them back to symbolic expressions.

I "cleaned up" your code, as you suggested.

You mentioned the possibility of using ithprime. Even better is nextprime. It's not really an optimization, but it leads to cleaner code.

Your idea of prime pairs with a fixed difference can be easily extended to prime constellations (aka tuples). For example, [11, 13, 17, 19] is a 4-tuple with difference pattern [2, 6, 8].

restart:
PrimeTuples:= proc(
    diffs::list(posint),
    #Three keyword-option stop criteria:
    {
        Interval::range({posint, infinity}):= 2..infinity,
        Timelimit::positive:= infinity,
        Count::{posint, infinity}:= infinity
    }
)
local
    P:= rtable(1..0), k:= 0, q:= lhs(Interval)-1,

    #Append next tuple to return array:
    Tuple:= proc()
    local T;
        q:= nextprime(q); T:= q +~ diffs;
        if andmap(isprime, T) then k:= k+1; P(k):= [q, T[]] fi;
        return
    end proc,

    #Check for cases with provably finitely many tuples: 
    Finite:= proc()
    local p:= 1, nd:= nops(diffs) + 1, df:= {diffs[]};
        do
            p:= nextprime(p);
            if p > nd then break fi;
            if nops(irem~(df, p) minus {0}) = p-1 then
                while q < p do Tuple() od; return true
            fi
        od;
        false
    end proc,

    #Check for cases with likely infinitely many tuples:        
    Infinite:= proc()
    option procname;
    local upper:= rhs(Interval);
        if {upper, Count, Timelimit} = {infinity} then
            error "no stop criterion"
        fi;
        while q < upper and k < Count do Tuple() od;
        return
    end proc       
;    
    if not Finite() then    
        try timelimit(eval(Timelimit, infinity= -1), Infinite())
        catch "time expired":
        end try
    fi;

    convert(P, list)
end proc
: 

PrimeTuples([2,4]);
                          [[3, 5, 7]]

PrimeTuples([2,8,14,26]);
            [[3, 5, 11, 17, 29], [5, 7, 13, 19, 31]]

PrimeTuples([2,4,6]);
                               []

PrimeTuples([2,6,8]);
Error, (in PrimeTuples) no stop criterion

PrimeTuples([2,6,8], Count= 9);
   [[5, 7, 11, 13], [11, 13, 17, 19], [101, 103, 107, 109], 
    [191, 193, 197, 199], [821, 823, 827, 829], [1481, 1483, 1487, 1489], 
    [1871, 1873, 1877, 1879], [2081, 2083, 2087, 2089], [3251, 3253, 3257, 3259]]

PrimeTuples([48], Interval= 2..1000);
[[5, 53], [11, 59], [13, 61], [19, 67], [23, 71], [31, 79], 
  [41, 89], [53, 101], [59, 107], [61, 109], [79, 127], 
  [83, 131], [89, 137], [101, 149], [103, 151], [109, 157], 
  [131, 179], [149, 197], [151, 199], [163, 211], [179, 227], 
  [181, 229], [191, 239], [193, 241], [223, 271], [229, 277], 
  [233, 281], [263, 311], [269, 317], [283, 331], [311, 359], 
  [331, 379], [349, 397], [353, 401], [373, 421], [383, 431], 
  [401, 449], [409, 457], [419, 467], [431, 479], [439, 487], 
  [443, 491], [461, 509], [499, 547], [509, 557], [521, 569], 
  [523, 571], [569, 617], [571, 619], [593, 641], [599, 647], 
  [613, 661], [643, 691], [653, 701], [661, 709], [691, 739], 
  [709, 757], [739, 787], [761, 809], [773, 821], [809, 857], 
  [811, 859], [829, 877], [839, 887], [859, 907], [863, 911], 
  [881, 929], [919, 967], [929, 977], [971, 1019], [983, 1031], 
  [991, 1039]] 

Download PrimeTuples.mw

There is a stock command for it: Fractals:-EscapeTime:-Newton. See its help page ?Fractals,EscapeTime,Newton.

The 2nd partial derivatives with respect to the 1st independent variable should be D[1,1](W)(0, tau) and D[1,1](W)(inf, tau).

The behavior you see is due to using an elementwise operator with multiple container arguments. It's not a bug, and has nothing to do with applyrule. To see what's happening, replace applyrule with foo.

The sublists of the list or Array of "neighbors" should not consist of vertex labels themselves but rather consist of the positive-integer indices into the vertex list that give those labels. So, this is correct:

GraphTheory:-Graph([0,1], [[2], [1]]);
      Graph 2: an undirected unweighted graph with 2 vertices and 1 edge(s)

I'll assume that that is sufficient explanation for you to correct the rest.

The only documentation of this is an extremely poorly written sentence in paragraph 6 of Description at ?GraphTheory,Graph:

• Note that the mapping is an integer mapping that indicates the vertices (if the vertices are labeled as  [1, 2, 3, ...]) or the location of each vertex in the vertex list V.

No example is given.

Also, the word "neighbors" is used incorrectly on this help page when referring to this argument to Graph; it should be "departures". Neighbors are departures plus arrivals. This distinction only matters for directed graphs.

Many graph theory problems can be solved more efficiently by working directly from the adjacency matrix. The following procedure is 20 times faster than your NumberOfK4.

restart:
CountCliques:= proc(G::GRAPHLN, q::posint)
uses GT= GraphTheory, It= Iterator, AT= ArrayTools;
local
    n:= GT:-NumberOfVertices(G),
    A:= AT:-Alias(
        1 -~ rtable(GT:-AdjacencyMatrix(G), datatype= integer[1]),
        [(0..n-1)$2]
    ),
    c, j, ct:= 0
;
    if q>2 then
        for j from 0 to n-1 do A[j,j]:= 0 od;
        for c in It:-Combination(n, q) do
            j:= [seq](c);
            if AT:-IsZero(A[j,j]) then ct++ fi
        od;
        ct
    elif q=1 then n
    else GT:-NumberOfEdges(G)
    fi
end proc
: 
GT:= GraphTheory:
G:= GT:-LineGraph(GT:-CompleteGraph(7)):
CountCliques(G, 3);
                              175

CodeTools:-Usage(CountCliques(G, 4));
memory used=5.73MiB, alloc change=0 bytes, 
cpu time=47.00ms, real time=49.00ms, gc time=0ns

                              105

You have several basic punctuation errors:

1. The last statement in the loop in step 8 ends with "HK." You need to remove the period.
2. You need to put a semicolon after "end do" in the following line.
3. At the end step 13, you need a semicolon after "end do".
4. In step 16, "NFLAG:= 0" should be "NFLAG = 0".
5. At the end of step 18, you need a semicolon after "end if".
6. In step 7, there is an assignment "W3:= W2 + HK*f (T,W2)". You have a blank space after f (it's very hard to see it in the 2D input). You need to remove the space.
7. Same problem with "f (T,W2)" in step 8.
8. Same problem with "f (T,W3)" in step 9.

At this point, after I make the above changes, I am unable to edit or execute the worksheet further. All GUI functions (within this particular worksheet) other than scrolling or saving are blocked; my other open worksheets are not affected. Saving and reopening the worksheet doesn't help.

I don't understand this at all, but I suspect that it has something to do with 2D Input. I've never seen this in the many tens of thousands of hours that I've used Maple over the past 20 years or so. I attach the worksheet here in case someone else can figure out what's going on. I don't think that this has anything to do with your coding; it's a GUI bug.

Download sol_odes.mw

It doesn't work because you don't have write permission to whatever directory (or folder) Maple considers to be your "current" directory. This is likely to be the main directory where the Maple system itself is stored, such as (on Windows) 

currentdir();
                 "C:\Program Files\Maple 2021"

Change that to something that you do have write permission for, such as your desktop. I (personally) would do

currentdir("C:/Users/carlj/desktop"): #Note **forward** slashes (/ not \)!

But you'd need to modify that because you're not "carlj". Once that's done, the save command in your worksheet should work.