This paper presents a methodology to study non-twist invariant circles and their bifurcations for area preserving maps, which is supported on the theoretical framework developed in Gonzalez-Enriquez et al. (Mem. Amer. Math. Soc. 227:vi+115, 2014). We recall that non-twist invariant circles are characterized not only by being invariant, but also by having some specified normal behavior. The normal behavior may endow them with extra stability properties (e.g., against external noise), and hence, they appear as design goals in some applications, e.g., in plasma physics, astrodynamics and oceanography. The methodology leads to efficient algorithms to compute and continue, with respect to parameters, non-twist invariant circles. The algorithms are quadratically convergent and, when implemented using FFT, have low storage requirement and low operations count per step. Furthermore, the algorithms are backed up by rigorous a posteriori theorems, proved and discussed in detail in Gonzalez-Enriquez et al. (Mem. Amer. Math. Soc. 227:vi+115, 2014), which give sufficient conditions guaranteeing the existence of a true non-twist invariant circle, provided an approximate invariant circle is known. Hence, one can compute confidently even very close to breakdown. With some extra effort, the calculations could be turned into computer-assisted proofs, see Figueras et al. (Found. Comput. Math. 17:1123–1193, 2017) for examples of the latter. The algorithms are also guaranteed to converge up to the breakdown of the invariant circles, and then, they are suitable to compute regions of parameters where the non-twist invariant circles exist. The calculations involved in the computation of the boundary of these regions are very robust, and they do not require symmetries and can run without continuous manual adjustments, largely improving methods based on the computation of very long period periodic orbits to approximate invariant circles. This paper contains a detailed description of our algorithms, the corresponding implementation and some numerical results, obtained by running the computer programs. In particular, we include calculations for two-dimensional parameter regions where non-twist invariant circles (with a prescribed frequency) exist. Indeed, we present systematic results in systems that do not contain symmetry lines, which seem to be unaccessible for previous methods. These numerical explorations lead to some open questions, also included here.