SUDOKU - rapidly find all possible solutions to a sudoku puzzle

Usage: Mout = sudoku(M)

M = intial sudoku matrix, with zeros for empty entries
Mout = solution as a 9x9 matrix if there is a unique solution, or as a 9x9xN matrix if there are N solutions

Notes:
(1) The algorithm employs recursion, but does as much as it can with straighforward deterministic deduction at each recursion level to increase overall speed.
(2) Supplying this function with an empty or overly sparse input matrix can create longer calulation times as the function searches for all possible solutions.
(3) A "No solution" error is generated if the input puzzle has no valid solution.
(4) Tested but no warranty, use at your own risk.
(5) Michael Kleder, December 2006

Examples:

% find the unique solution to this puzzle in a fraction of a second:
M = [0 0 8 0 9 0 5 0 0;0 0 1 0 7 0 4 0 0;0 0 4 0 3 0 6 0 0;
0 1 0 0 0 6 0 0 7;0 9 0 0 0 3 0 0 0;0 2 0 0 5 0 0 6 0;
0 5 0 0 4 0 0 2 0;0 0 0 8 0 0 0 3 0;6 0 0 1 0 0 0 4 0];
sudoku(M)

% find the 100 possible solutions to this puzzle in few seconds:
M = [0 0 8 0 9 0 5 0 0;0 0 1 0 7 0 0 0 0;0 0 4 0 3 0 6 0 0;
0 1 0 0 0 0 0 0 7;0 9 0 0 0 3 0 0 0;0 2 0 0 5 0 0 6 0;
0 5 0 0 4 0 0 2 0;0 0 0 8 0 0 0 3 0;0 0 0 1 0 0 0 4 0];
tic;M=sudoku(M);toc;size(M,3)