Abstract. In the paper, the method of penalized high order Legendre polynomials and specific designs of experiments for model building on the basis of numerical or physical experiments is proposed. The use of multivariate orthogonal Legendre polynomials allows relatively simple elimination of insignificant terms. In this paper the generalized thin-plate energy functional is used for penalization of the least square functional. The optimal choice of smoothing parameter (penalization coefficient) is implemented using the cross-validation method. Special axial-symmetrical and 90 degree rotational symmetrical D-optimal experimental designs are obtained using the direct constrained optimization method. The proposed algorithm is implemented in the software tool EDAOpt that has been developed in the Riga Technical University for more than 13 years [1, 2]. The method was tested for known optimization test problems with 2-5 variables and showed prediction accuracy comparable with kriging. For noisy responses the proposed approach gives better prediction accuracy of the approximate model.