Sometimes---more often than I would expect---solving a more general problem is easier yet still gives a complete solution to the original problem. Such is the case here. The problem at hand is equivalent to finding all directed paths of a given length in a graph. What you call routes are formally called paths: a sequence of distinct vertices, each adjacent to the preceding.

Once the data structure is defined as I have done in the procedure below, finding the paths is essentially a one-liner. There's no need to consider special cases like vertices on corners, sides, etc.

**AllPaths:= proc(**

** G::GraphTheory:-Graph**

** ,maxlength::nonnegint:= GraphTheory:-NumberOfVertices(G)**

**)**

**#Returns a table R such that R[k] is a list of all directed paths of length k in G,**

**#where k ranges from 0 to maxlength. G can be directed or undirected.**

**#Each path is represented as a list of vertices.**

**#The algorithm is breadth-first search.**

** uses GT= GraphTheory;**

** local**

** k #current length **

** ,p #current path**

** ,v #current vertex**

** ,V:= GT:-Vertices(G) #constant**

** ,E:= table(V =~ (`{}`@op)~ (GT:-Departures(G))) #constant table of edges**

** ,R:= table([0= `[]`~ (V)]) #Initialize with zero-length paths.**

** ;**

** for k to maxlength do**

** #For each path of length k-1, make all possible 1-vx extensions.**

** R[k]:= [seq](seq([p[], v], v in E[p[-1]] minus {p[]}), p in R[k-1])**

** end do; **

** eval(R) **

**end proc:**

Now apply it to the problem at hand:

**squares:= [seq](seq([i,j], i= 1..4), j= 1..4):**

Name the vertices conveniently.

**V:= (v-> cat(['a','b','c','d'][v[1]], v[2]))~ (squares);**

V := [a1, b1, c1, d1, a2, b2, c2, d2, a3, b3, c3, d3, a4, b4, c4, d4]

**Vsq:= table(V =~ squares):**

Define the edges of the graph. Horizontal, vertical, and diagonal adjacency are all covered by simply saying that the infinity norm of the difference of the points is 1.

**Sq44:= GraphTheory:-Graph(**

select(e-> LinearAlgebra:-Norm(<Vsq[e[1]] - Vsq[e[2]]>) = 1, combinat:-choose({V[]}, 2))

):

Generate the paths and time it.

**CodeTools:-Usage(assign(Pths, AllPaths(Sq44, 6))):**

memory used=11.31MiB, alloc change=4.27MiB, cpu time=219.00ms, real time=210.00ms

Explore the results a bit.

**nops(Pths[6]);**

68272

Get a well-distributed sample of the paths.

**Samp:= [seq](Pths[6][2^k], k= 1..16);**

Samp := [[a1, a2, a3, a4, b3, b2, c1], [a1, a2, a3, a4, b3, b2, c3], [a1, a2, a3, a4, b3, c2, b2], [a1, a2, a3, a4, b3, c3, c2], [a1, a2, a3, a4, b4, c3, c4], [a1, a2, a3, b2, b3, c3, d4], [a1, a2, a3, b2, c3, d3, c4], [a1, a2, a3, b4, c3, b3, c4], [a1, a2, b1, c2, c3, b4, a3], [a1, a2, b3, b4, a3, b2, c2], [a1, b1, c1, c2, d3, c3, b3], [a2, a1, b2, a3, b3, c2, d1], [a2, b3, c3, d2, c1, c2, b2], [a4, b4, c3, c2, b1, c1, b2], [b4, c3, d2, c1, b1, a2, a1], [d4, c3, d2, c2, b3, b4, a3]]

Develop a few tools for visualizing the paths.

To center a vertex in its square:

**pt:= v-> Vsq[v] -~ 1/2:**

Background to plot on:

**grid:= plot(0, view= [0..4, 0..4], gridlines, scaling= constrained):**

**pathplot:= p->**

plots:-display(

grid

,plot(

[[pt](p[1]), [pt](p[-1]), pt~(p)]

,style= [point, point, line], color= [green, red, yellow]

,symbol= diamond, symbolsize= 30

)

)

:

And finally here are the 16 sample paths:

**plots:-display(Matrix(4,4, pathplot~ (Samp)));**

(Sorry, I don't know how to upload an array of plots. If you execute the above code yourself, you will get a 4x4 array of plots---one for each of the sample paths.)