  Using highly efficient GPU-based simulations of the tight-binding
Bogoliubov-de Gennes equations we solve self-consistently for the pair
correlation in rhombohedral (ABC) and Bernal (ABA) multilayer graphene by
considering a finite intrinsic s-wave pairing potential. We find that the two
different stacking configurations have opposite bulk/surface behavior for the
order parameter. Surface superconductivity is robust for ABC stacked multilayer
graphene even at very low pairing potentials for which the bulk order parameter
vanishes, in agreement with a recent analytical approach. In contrast, for
Bernal stacked multilayer graphene, we find that the order parameter is always
suppressed at the surface and that there exists a critical value for the
pairing potential below which no superconducting order is achieved. We
considered different doping scenarios and find that homogeneous doping strongly
suppresses surface superconductivity while non-homogeneous field-induced doping
has a much weaker effect on the superconducting order parameter. For multilayer
structures with hybrid stacking (ABC and ABA) we find that when the thickness
of each region is small (few layers), high-temperature surface
superconductivity survives throughout the bulk due to the proximity effect
between ABC/ABA interfaces where the order parameter is enhanced.
