> 
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:

> 
Numbrix := proc( M :: ~Matrix, { inline :: truefalse := false } )
local S, adjacent, eq, i, initial, j, k, kk, m, n, one, single, sol, unique, val, var, x;
(m,n) := upperbound(M);
initial := &and(seq(seq(ifelse(M[i,j] = 0
, NULL
, x[i,j,M[i,j]]
)
, i = 1..m)
, j = 1..n));
adjacent := &and(seq(seq(seq(x[i,j,k] &implies &or(NULL
, ifelse(i>1, x[i1, j, k+1], NULL)
, ifelse(i<m, x[i+1, j, k+1], NULL)
, ifelse(j>1, x[i, j1, k+1], NULL)
, ifelse(j<n, x[i, j+1, k+1], NULL)
)
, i = 1..m)
, j = 1..n)
, k = 1 .. m*n1));
one := &or(seq(seq(x[i,j,1], i=1..m), j=1..n));
single := ¬(&or(seq(seq(seq(seq(x[i,j,k] &and x[i,j,kk], kk = k+1..m*n), k = 1..m*n1)
, i = 1..m), j = 1..n)));
sol := Logic:Satisfy(&and(initial, adjacent, one, single));
if sol = NULL then
error "no solution";
end if;
if inline then
S := M;
else
S := Matrix(m,n);
end if;
for eq in sol do
(var, val) := op(eq);
if val then
S[op(1..2, var)] := op(3,var);
end if;
end do;
S;
end proc:

Two simple examples
> 
A:=<0,0,5; 0,0,0; 0,0,9>;
# The unique solution
NumbrixPuzzle(A);
A:=<0,0,5; 0,0,0; 0,8,0>;
# 4 solutions
NumbrixPuzzle(A);


(1) 
Comparison with Numbrix procedure. The example is taken from
http://rosettacode.org/wiki/Solve_a_Numbrix_puzzle
> 
A:=<0, 0, 0, 0, 0, 0, 0, 0, 0;
0, 0, 46, 45, 0, 55, 74, 0, 0;
0, 38, 0, 0, 43, 0, 0, 78, 0;
0, 35, 0, 0, 0, 0, 0, 71, 0;
0, 0, 33, 0, 0, 0, 59, 0, 0;
0, 17, 0, 0, 0, 0, 0, 67, 0;
0, 18, 0, 0, 11, 0, 0, 64, 0;
0, 0, 24, 21, 0, 1, 2, 0, 0;
0, 0, 0, 0, 0, 0, 0, 0, 0>;
CodeTools:Usage(NumbrixPuzzle(A));
CodeTools:Usage(Numbrix(A));

memory used=7.85MiB, alloc change=3.01MiB, cpu time=172.00ms, real time=212.00ms, gc time=93.75ms


memory used=1.21GiB, alloc change=307.02MiB, cpu time=37.00s, real time=31.88s, gc time=9.30s



(2) 
In the example below, which has 104 solutions, the Numbrix procedure is faster.
> 
C:=Matrix(5,{(1,1)=1,(5,5)=25});
CodeTools:Usage(NumbrixPuzzle(C)):
nops([%]);
CodeTools:Usage(Numbrix(C)):

memory used=0.94GiB, alloc change=22.96MiB, cpu time=12.72s, real time=11.42s, gc time=2.28s


memory used=34.74MiB, alloc change=0 bytes, cpu time=781.00ms, real time=783.00ms, gc time=0ns


