In this paper, a new numerical procedure based on a cycle-by-cycle analysis has been constructed for creep-fatigue behavior and life prediction of high-temperature structures under multi-axial stress states. Within this framework, a modified unified viscoplastic constitutive model with isotropic hardening and modified kinematic hardening rules is developed to simulate the cycle-by-cycle stress-strain responses. Moreover, the newly constructed creep-fatigue approach calculates fatigue and creep damage variables using the critical plane method (CPM) and the modified strain energy density exhaustion (SEDE) model, respectively. The multi-axial ductility factor and elastic follow-up factor are also introduced into the modified SEDE model to accommodate the special multi-axial and mixed controlled modes, which are widely existed in practical structures. In order to validate the feasibility of the proposed numerical procedure, a series of creep-fatigue tests of notched specimens made from nickel-based GH4169 superalloy were carried out at 650 °C. The predicted numbers of cycles to crack initiation agree well with the experimental data. Evidence of crack initiation under various loading conditions was observed via the electron backscatter diffraction (EBSD) technique, indicating location-dependent crack initiations depending on loading conditions. In detail, the crack initiation sites shifting from surface to subsurface with increasing hold times can be well simulated by the proposed numerical procedure due to a reasonable description of the creep-fatigue damage evolution.