PSR J2021+4026 in the Cygnus region is one of the most interesting radio-quiet pulsars detected by Fermi-LAT. This source, thought to be associated with the Gamma Cygni supernova remnant, is the first isolated gamma-ray pulsar that exhibited clear evidence of variability, with a simultaneous flux and spindown change first occurred in October 2011. After a long recovery phase, in August 2015 the pulsar reached the pre-2011 flux and timing characteristics and in September 2017 underwent a new mode change. This behaviour is still unique among the population of gamma-ray pulsars and therefore we performed multiwavelength follow-up observations in order to understand the physics behind this event. We have analyzed Fermi-LAT gamma rays simultaneous to two deep XMM-Newton observations taken after the flux drop and during the recovery phase. We also analyzed a deep observation with the Gran Telescopio Canarias carried on after the recovery in order to search for its optical counterpart. We present the results of this multiwavelength campaign and how such studies can help constraining the pulsar geometry and deciphering these rare gamma-ray variability events.