> 
SerpentinePaths:=proc(S::list, P::list:=[1,1])
local OneStep, A, m, F, B, T, a;
OneStep:=proc(A::listlist)
local s, L, B, T, k, l;
s:=max[index](A);
L:=[[s[1]1,s[2]],[s[1]+1,s[2]],[s[1],s[2]1],[s[1],s[2]+1]];
T:=table(); k:=0;
for l in L do
if l[1]>=1 and l[1]<=S[1] and l[2]>=1 and l[2]<=S[2] and A[op(l)]=0 then k:=k+1; B:=subsop(l=a+1,A);
T[k]:=B fi;
od;
convert(T, list);
end proc;
A:=convert(Matrix(S[], {(P[])=1}), listlist);
m:=S[1]*S[2]1;
B:=[$ 1..m];
F:=LM>ListTools:FlattenOnce(map(OneStep, `if`(nops(LM)<=30000,LM,LM[30000..1])));
T:=[A];
for a in B do
T:=F(T);
od;
map(convert, T, Matrix);
end proc:

> 
NumbrixPuzzle:=proc(A::Matrix)
local A1, L, N, S, MS, OneStepLeft, OneStepRight, F1, F2, m, L1, p, q, a, b, T, k, s1, s, H, n, L2, i, j, i1, j1, R;
uses ListTools;
S:=upperbound(A); N:=nops(op(A)[3]); MS:=`*`(S);
A1:=convert(A, listlist);
for i from 1 to S[1] do
for j from 1 to S[2] do
for i1 from i to S[1] do
for j1 from 1 to S[2] do
if A1[i,j]<>0 and A1[i1,j1]<>0 and abs(A1[i,j]A1[i1,j1])<abs(ii1)+abs(jj1) then return `no solutions` fi;
od; od; od; od;
L:=sort(select(e>e<>0, Flatten(A1)));
L1:=[`if`(L[1]>1,seq(L[1]k, k=0..L[1]2),NULL)];
L2:=[seq(seq(`if`(L[i+1]L[i]>1,L[i]+k,NULL),k=0..L[i+1]L[i]2), i=1..nops(L)1), `if`(L[1]<MS,seq(L[1]+k,k=0..MSL[1]1),NULL)];
OneStepLeft:=proc(A1::listlist)
local s, M, m, k, T;
uses ListTools;
s:=Search(a, Matrix(A1));
M:=[[s[1]1,s[2]],[s[1]+1,s[2]],[s[1],s[2]1],[s[1],s[2]+1]];
T:=table(); k:=0;
for m in M do
if m[1]>=1 and m[1]<=S[1] and m[2]>=1 and m[2]<=S[2] and A1[op(m)]=0 then k:=k+1; T[k]:=subsop(m=a1,A1);
fi;
od;
convert(T, list);
end proc;
OneStepRight:=proc(A1::listlist)
local s, M, m, k, T, s1;
uses ListTools;
s:=Search(a, Matrix(A1)); s1:=Search(a+2, Matrix(A1));
M:=[[s[1]1,s[2]],[s[1]+1,s[2]],[s[1],s[2]1],[s[1],s[2]+1]];
T:=table(); k:=0;
for m in M do
if m[1]>=1 and m[1]<=S[1] and m[2]>=1 and m[2]<=S[2] and A1[op(m)]=0 and `if`(a+2 in L, `if`(is(abs(s1[1]m[1])+abs(s1[2]m[2])>1),false,true),true) then k:=k+1; T[k]:=subsop(m=a+1,A1);
fi;
od;
convert(T, list);
end proc;
F1:=LM>ListTools:FlattenOnce(map(OneStepLeft, LM));
F2:=LM>ListTools:FlattenOnce(map(OneStepRight, LM));
T:=[A1];
for a in L1 do
T:=F1(T);
od;
for a in L2 do
T:=F2(T);
od;
R:=map(t>convert(t,Matrix), T);
if nops(R)=0 then return `no solutions` else R fi;
end proc:

Simple examples
> 
SerpentinePaths([3,3]); # All the serpentine paths for the matrix 3x3, starting with [1,1]position
SerpentinePaths([3,3],[1,2]); # No solutions if the start with [1,2]position
SerpentinePaths([4,4]): # All the serpentine paths for the matrix 4x4, starting with [1,1]position
nops(%);
nops(SerpentinePaths([4,4],[1,2])); # The number of all the serpentine paths for the matrix 4x4, starting with [1,2]position
nops(SerpentinePaths([4,4],[2,2])); # The number of all the serpentine paths for the matrix 4x4, starting with [2,2]position


(1) 
Below we find 12,440 serpentine paths in the matrix 6x6 starting from various positions (the set L )
> 
k:=0: n:=6:
for i from 1 to n do
for j from i to n do
k:=k+1; S[k]:=SerpentinePaths([n,n],[i,j])[];
od: od:
L1:={seq(S[i][], i=1..k)}:
L2:=map(A>A^%T, L1):
L:=L1 union L2:
nops(L);


(2) 
Further, using the list L, we generate 20 examples of Numbrix puzzles with the unique solutions
> 
T:='T':
N:=20:
M:=[seq(L[i], i=combinat:randcomb(nops(L),N))]:
for i from 1 to N do
for k from floor(n^2/4) do
T[i]:=Matrix(n,{seq(op(M[i])[3][j], j=combinat:randcomb(n^2,k))});
if nops(NumbrixPuzzle(T[i]))=1 then break; fi;
od: od:
T:=convert(T,list):
T1:=[seq([seq(T[i+j],i=1..5)],j=[0,5,10,15])]:
DocumentTools:Tabulate(Matrix(4,5, (i,j)>T1[i,j]), fillcolor = "LightYellow", width=95):

The solutions of these puzzles
> 
DocumentTools:Tabulate(Matrix(4,5, (i,j)>NumbrixPuzzle(T1[i,j])[]), fillcolor = "LightYellow", width=95):

