The nonlinear three-dimensional magnetohydrodynamics (MHD) code NFTC and the basic algorithms included in the code are presented. NFTC code is used for simulations of nonlinear neoclassical instabilities (or neoclassical tearing modes (NTMs)) in tokamaks. An effective fully implicit numerical scheme allows the transport profile to evolve self-consistently with the nonlinear MHD instabilities and externally applied sources and sinks. It is shown that numerical code has enough accuracy for simulation of plasma instabilities in arbitrary three-dimensional geometry. The coefficients of a special finite-difference scheme in the radial direction are constructed in the form which gives the exact asymptotic dependence of the solution by the ratio of small parameters near high derivatives and the mesh step size. The reproduction of the nonlinear threshold effects for NTMs in conditions of anisotropic heat conductivity is presented.