Solve a magic square variant in Matlab

fu_puzzle.m

% Program to solve a sort of magic square, where only the rows and columns
% must sum to the magic constant, and not the diagonals.  Some elements of
% the square are given in advance while others (indicated with 0) are
% unknown.
%
% TF Time-Wasting Project 04-2013

function result = fu_puzzle
puzzle = ...
  [ 1  0  0  0 22
    0 24  0 11  0
    0  0  5  0  0
    0  3  0 17  0
   19  0  0  0 13 ];
result = solve_magic(puzzle);
end

function M = magic_constant(n)
M = n*(n^2 + 1)/2;
end

function solvable = is_solvable_helper(puzzle, n, M)

solvable = 0;

% solvable = all(   ((~any(puzzle==0)) & (sum(puzzle) == M)) ...
%   | (( any(puzzle==0)) & (sum(puzzle) <=  M)) );

values = remaining_values(puzzle);
for ii=1:n
  row = puzzle(ii, :);   % Fetch a row of the puzzle
  x = (row==0);          % Identify the unknown cells
  nx = numel(find(x));   % How many are they?
  
  % If all cells in the row are filled in, then the row MUST sum to the
  % correct value.
  if (nx == 0) && (sum(row) ~= M) 
    return
  end
  
  % If not all cells are filled in, populate them with the lowest-valued
  % remaining possibilities, and ensure that the sum of the row does not
  % exceed the magic constant.
  row(x) = values(1:nx);
  if sum(row) > M
    return
  end
  
  row(x) = values(end-nx+1:end);
  if sum(row) < M
      return
  end
end
% If we get this far, then we haven't found any problems.
solvable = 1;
end


% Check whether a puzzle may still be solved by filling in more values
function result = is_solvable(puzzle)
n = size(puzzle, 1);
M = magic_constant(n);
result = is_solvable_helper(puzzle, n, M) && ...
         is_solvable_helper(puzzle', n, M);
end

% What numbers haven't been used yet?
function values = remaining_values(puzzle)
values = setdiff(1:numel(puzzle), puzzle(:));
end

% Solve the puzzle
function solutions = solve_magic(puzzle)

solutions = {};
if ~is_solvable(puzzle)  
  return
end

values = remaining_values(puzzle);

if isempty(values)
  % No more values to put in --> the problem is solved!
  solutions = {puzzle}
  return
end

%position = find(puzzle==0, 1, 'first');
 
% Try to encounter a constraint as soon as possible:
positions = find(puzzle==0);
[rows, cols] = ind2sub(size(puzzle), positions);
col_remaining = sum(puzzle ==0);
row_remaining = sum(puzzle'==0);
remaining = sort([row_remaining(rows); col_remaining(cols)]);
[~, ii] = sortrows(remaining');
position = positions(ii(1));

for value = values,
  puzzle(position) = value;
  solutions = [solutions solve_magic(puzzle)];
end
end