Simultaneous observations of PSR B0823+26 with ESAs XMM-Newton, the Giant Metrewave Radio Telescope and international stations of the Low Frequency Array revealed synchronous X-ray/radio switching between a radio-bright (B) mode and a radio-quiet (Q) mode. During the B mode we detected PSR B0823+26 in 0.2$-$2 keV X-rays and discovered pulsed emission with a broad sinusoidal pulse, lagging the radio main pulse by 0.208 $pm$ 0.012 in phase, with high pulsed fraction of 70$-$80%. During the Q mode PSR B0823+26 was not detected in X-rays (2 $sigma$ upper limit a factor ~9 below the B-mode flux). The total X-ray spectrum, pulse profile and pulsed fraction can globally be reproduced with a magnetized partially ionized hydrogen atmosphere model with three emission components: a primary small hot spot ($T$$sim$3.6$times10^6$ K, $R$$sim$17 m), a larger cooler concentric ring ($T$$sim$1.1$times10^6$ K, $R$$sim$280 m) and an antipodal hot spot ($T$$sim$1.1$times10^6 $ K, $R$$sim$100 m), for the angle between the rotation axis and line of sight direction $sim66^circ$. The latter is in conflict with the radio derived value of $(84pm0.7)^circ$. The average X-ray flux within hours-long B-mode intervals varied by a factor $pm$20%, possibly correlated with variations in the frequency and lengths of short radio nulls or short durations of weak emission. The correlated X-ray/radio moding of PSR B0823+26 is compared with the anti-correlated moding of PSR B0943+10, and the lack of X-ray moding of PSR B1822-09. We speculate that the X-ray/radio switches of PSR B0823+26 are due to variations in the rate of accretion of material from the interstellar medium through which it is passing.