The aim of the present paper is to present numerical procedures, using two types of FE methods, for the secondary and higher order bifurcations (tertiary and quarternary bifurcations) which appear after the first buckling of thin shell structures. With a limitation to a cylindrical shell under axial compression, the post-bifurcation behavior after the first bifurcation is analyzed precisely by considering modal coupling between several deformation modes of higher order harmonic wave numbers, and on all the way of post-bifurcation path the positive definiteness of incremental stiffness matrix of uncoupled modes is examined step by step to find if any higher order bifurcation or transcritical bifurcation appears. Based on the step-by-step examination, the higher order bifurcation points and the corresponding bifurcation modes are traced.