In your script file, you did a spin-polarized calculation for the ferromagnetic state of zigzag graphene nanoribbon. In addition, 4 x 1 supercell of zigzag graphene nanoribbon is used. So, the band dispersion of two bands near the Fermi levels are folded in your calculation results. This makes your results different from the published one, which is obtained by the primitive unit cell of zigzag graphene nanoribbon.
Please do the anti-ferromagnetic state of zigzag graphene nanoribbon using its primitive unit cell. You will reproduce the published result exactly.