We describe a new multifractal finite size scaling (MFSS) procedure and its application to the Anderson localization-delocalization transition. MFSS permits the simultaneous estimation of the critical parameters and the multifractal exponents. Simulations of system sizes up to L^3=120^3 and involving nearly 10^6 independent wavefunctions have yielded unprecedented precision for the critical disorder W_c=16.530 (16.524,16.536) and the critical exponent nu=1.590 (1.579,1.602). We find that the multifractal exponents Delta_q exhibit a previously predicted symmetry relation and we confirm the non-parabolic nature of their spectrum. We explain in detail the MFSS procedure first introduced in our Letter [Phys. Rev. Lett. 105, 046403 (2010)] and, in addition, we show how to take account of correlations in the simulation data. The MFSS procedure is applicable to any continuous phase transition exhibiting multifractal fluctuations in the vicinity of the critical point.