Codes based on SUDOKU puzzles are discussed, and belief propagation decoding introduced for the erasure channel. Despite the non-linearity of the code constraints, it is argued that density evolution can be used to analyse code performance due to the invariance of the code under alphabet permutation. The belief propagation decoder for erasure channels operates by exchanging messages containing sets of possible values. Accordingly, density evolution tracks the probability mass functions of the set cardinalities. The equations governing the mapping of those probability mass functions are derived and calculated for variable and constraint nodes, and decoding thresholds are computed for long SUDOKU codes with random interleavers.