You have mixed up the monolayer and bilayer input files, first of all. Second, your bilayer structure is not really a bilayer structure but rather an infinite stack of MoS2 layers (so, actually it's bulk MoS2). You would need to make more vacuum in the C direction to make it bilayer. Also, there is no reason to have k-point sampling in the C direction, it only makes the calculations take more time (your "monolayer" file which is really not a monolayer has this setting).
Do you have a reference for the picture from the article? Otherwise it's hard to compare, one needs to look into the precise details also of their calculation. I have seen comments about LDA/GGA being insufficient to treat some long-range interactions in MoS2, but I haven't had time to look into the details.