The spin system of the 2d Ising model having a hexagonal-lattice is simulated using non-deterministic Cellular Automata (CA). The method to implement this program is outlined and our results show a good approximation to the exact analytic solution for the initial random lattice. The 2d system is studied with a 40 × 40 hexagonal-lattice with five different boundary conditions (bcs) i.e., adiabatic, periodic, reflexive, fixed +1 and fixed -1 with random orientation of spins as initial conditions in the absence of an external applied magnetic field. The critical temperature below which the spontaneous magnetization appears as well as other physical quantities such as the magnetisation, energy, specific heat, susceptibility and entropy with each of the bcs are calculated. The phase transition occurs around T H c = 1.5 which approximates well with the result obtained from exact analytic solution by Wannier and Houtappel. The computational advantage of using non-deterministic Cellular Automata is discussed.