Quantum algorithms are good at solving combinatorial situations, a challenge faced by many NP-hard and NP-complete optimization problems. In general, each non-polynomial (NP) problem has its own characteristics and requires special attention. Despite this fact, one hard problem can naturally be deduced from others and the approaches of solving one may be helpful to solve others. Therefore, this study develops a three-dimensional quantum Ising model to solve the NP-Complete Sudoku puzzle games, based on the widely adopted approach of optimizing Hamiltonian energy to model NP-hard Ising spin glasses. This developed paradigm can benefit software practitioners to solve more NP problems. In our proposed methodology, mathematical foundations are first formulated based on the couplings of atomic spins to model the puzzle rules and constraints. Next, quantum modeling is discussed to demonstrate the construction of observable operators from Pauli gates to compute the expectation values of plausible quantum states for the puzzle. Finally, individual quantum algorithms are integrated to become a Sudoku solver. The solver is investigated through the execution of quantum approximate optimization algorithm (QAOA) combined with the constrained optimization by linear approximation (COBYLA) optimizer provided in the IBM Qiskit SDK. A code snippet for one of the puzzle rules is presented to verify the outcome of the developed quantum circuits.