Aiming at the problems of identification difficulties and low identification accuracy in modelling and identification of multiple-input multiple-output (MIMO) nonlinear Gaussian time-varying systems, this paper proposes an identification scheme based on the step-by-step approximation of multidimensional Taylor network (MTN). The aim of this paper is to improve the modelling of complex nonlinear systems so as to improve the prediction performance and control effect of the system. Different from the traditional multidimensional Taylor network identification method, this method adopts an order-by-order approximation strategy, which seeks its parameters sequentially from the lower order to the higher order, and continuously optimises the parameter weights during the parameter seeking process. Firstly, the nonlinear function model is approximated as a polynomial form by the order-by-order Taylor expansion, and then the weight parameters of each order of the Taylor expansion are calculated and updated step by step by using the algorithm based on the Variable Forgetting Factor Recursive Least Squares (VFF-RLS) method. Through iterative optimized of these parameters, dynamic weight assignment to each order of the Taylor expansion is achieved. A parameter-identified nonlinear function model is finally obtained, which can more accurately describe the dynamic behaviour and characteristics of the system. Finally, an arithmetic simulation is carried out through an example to verify the effectiveness of the proposed method.