It’s easiest to check the finite differences for second derivatives against the Hessian (which is what I do in GULP when implementing models) as there are no phase factors to worry about.
For graphene be careful about comparing eigenvectors - if you have degenerate modes (which is highly likely in high symmetry cases) then the answers can vary since the modes can be mixed anyway you like & so there is not necessarily a unique solution. Always test for low symmetry cases to avoid this.