f (R) theory is a modification of Einstein general relativity which has many interesting results in cosmology and astrophysics. To derive solutions for black holes solution in this theory is difficult due to the fact that it is fourth order differential equations. In this study, we use the quadratic form of f (R) = R + βR 2 − 2Λ, where β is a dimensional parameter and Λ is the cosmological constant, to find the solutions of new 4-dimension spherically symmetric non-charged and charged NUT and Taub-NUT black holes. An interesting property for these black holes is the fact that NUT spacetime has a non-vanishing value of the magnetic field, while Taub-NUT one does magnetic and electric fields. We calculate the energy conditions of the charged black holes and show that all of them are satisfied for the Taub-NUT spacetime. Finally, we study some thermodynamic quantities such as entropy, temperature, specific heat and Gibbs free energy. The calculations of heat capacity and free energy show that the charged NUT and Taub NUT black holes have local stability.Recently, observations confirm that our universe suffers from acceleration. Since that, the correspondence between the accelerated epoch and the inflationary mechanism help scientists to assume that dark energy may have a geometric origin [3][4][5]. Later, scientists discovered that the terms which are related to the quantum higher order corrections are responsible for the accelerating expansion rate of our universe at large structures. This phenomenon has been investigated at the fourth order of f (R) [6][7][8][9][10][11][12][13] (for more references on modified gravity theories as well as dark energy problem, see, e.g., [14,15]).The mechanism which determines the difference between f (R) gravitational theories and the Einstein general relativity (GR), that is ensured to be detected at the scales of astrophysics and strong gravitational field, is to derive black holes that are different from those of GR (vacuum or electro-vacuum) [16][17][18][19][20][21]. There are many applications carried out on f (R), among them are: The one-loop effective action of f (R) theories on the de-Sitter background that has been explained in [22]. One of the most interesting phenomena which can occur in charged black hole solutions is the anti-evaporation process which explains the primordial black hole. In order that the anti-evaporation may occur in the Einstein theory, one needs to involve the quantum correction terms from the matter field. Anti-evaporation could occur at the classical level on modified gravity theories like f (R) gravity [23,24]. A subject to acquire solutions for static spherically symmetric black holes in f (R) gravity has been addressed in [25]. In f (R) gravity, solutions for spherically symmetric black holes have been obtained by assuming that the scalar curvature is constant [26][27][28]. Also, spherically symmetric non-charged and charged black hole solutions have been discussed with no constraints on the scalar curvature and the Ricci tensor in [29][30][31]....