A multi-scale MHD numerical scheme is developed for analysis of nonlinear evolution of a beta-increasing plasma. The scheme is based on iterative calculations of nonlinear dynamics based on the reduced MHD (RMHD) equations and three-dimensional static equilibrium. The equation for average pressure in the RMHD equations plays the role of a transport equation that involves a heat source term and background pressure diffusion terms. The heat source term is controlled so that the beta value should be increased at a constant rate. The scheme is applied to a Large Helical Device (LHD) plasma up to average beta of 1.05%, which is unstable against linear ideal interchange modes while beta values much higher than the stability limit are obtained in the experiments. The result with the multi-scale scheme indicates that many local flat regions are generated in the background pressure profile in the nonlinear evolution of the interchange modes. This structure of the pressure profile suppresses disruptive phenomena because it reduces the driving force of the modes at higher beta value. Such self-organization in the pressure profile is considered to be the stabilizing mechanism in the plasma.