Variational and numerical study of solution switching in two and three core couplers with higher order perturbations