In statistical analyses, especially those using a multiresponse regression model approach, a mathematical model that describes a functional relationship between more than one response variables and one or more predictor variables is often involved. The relationship between these variables is expressed by a regression function. In the multiresponse nonparametric regression (MNR) model that is part of the multiresponse regression model, estimating the regression function becomes the main problem, as there is a correlation between the responses such that it is necessary to include a symmetric weight matrix into a penalized weighted least square (PWLS) optimization during the estimation process. This is, of course, very complicated mathematically. In this study, to estimate the regression function of the MNR model, we developed a PWLS optimization method for the MNR model proposed by a previous researcher, and used a reproducing kernel Hilbert space (RKHS) approach based on a smoothing spline to obtain the solution to the developed PWLS optimization. Additionally, we determined the symmetric weight matrix and optimal smoothing parameter, and investigated the consistency of the regression function estimator. We provide an illustration of the effects of the smoothing parameters for the estimation results using simulation data. In the future, the theory generated from this study can be developed within the scope of statistical inference, especially for the purpose of testing hypotheses involving multiresponse nonparametric regression models and multiresponse semiparametric regression models, and can be used to estimate the nonparametric component of a multiresponse semiparametric regression model used to model Indonesian toddlers' standard growth charts.