Improved double-integration technique to approximate integral-balance solutions of non-linear fractional subdiffusion equations has been conceived. The time-fraction subdiffusion equation with Dirichlet boundary condition and a powerlaw fractional diffusivity has been chosen as a test example. Problems pertinent to approximation of time-fractional Riemann-Liouville derivative when the distribution is expressed as a parabolic profile with unspecified exponent and accuracy of the solutions have been analyzed. The final solution is a closed-form can be presented with either a similarity variable of a fractional order as independent variable or by an effective similarity variable incorporating the effects of both the fractional order and the nonlinearity of the diffusion coefficient. Optimization problem pertinent to determination of the optimal exponent of the parabolic profile, dependent on both the fractional order and the nonlinearity parameter of the diffusion coefficient, has been developed by a modified technique transforming the time-varying domain of integration into one with fixed boundaries. It was clearly defined that the approximate profile can exhibit a concave behaviour, typical for subdiffusion relaxation processes when the non-linearity of the diffusion coefficient is low and the fractional order is high. Otherwise, the increase in the nonlinearity of the diffusion coefficient results in convex profiles typical for the degenerate diffusion behaviour.