This study presents a comprehensive framework for fitting the phenomenological potential model parameters of charmonium by integrating neural network wavefunctions. This approach overcomes the limitations of traditional fitting methods, which adapt slowly to new experimental data. In this method, the radial Schrödinger equation is solved variationally using a neural network wavefunction ansatz. The wavefunction incorporates an explicitly constructed node factor that accurately represents the nodal structure of excited states. It also employs a residual network architecture to enhance expressiveness and training stability. The number of nodes is specified according to the target radial excitation, while the node positions are treated as trainable parameters, enabling the neural network to optimize the excited-state wavefunctions without preassigning their nodal locations. The energy variance is used as the training objective, and automatic differentiation is employed to calculate the derivatives required in the variational optimization and subsequent parameter inference. After determining the eigenwavefunctions and eigenenergies, the potential parameters are inferred via maximum likelihood estimation with the Laplace approximation, yielding optimal parameter values and their uncertainties. In the likelihood construction, a fixed theoretical uncertainty is introduced together with the experimental errors so that the fitting procedure reflects both experimental precision and the intrinsic approximation of the potential model. The method is applied to a phenomenological charmonium potential that includes a Coulomb term, a linear confinement term, a smeared spin-spin interaction, and spin-orbit and tensor forces treated perturbatively. The fitted parameters are the charm quark mass _mc = m_c = 1.4801 ±0.0055 GeV, the strong coupling constant α_s = 0.5459 ±0.0120, the string tension σ_lin = 0.1442 ±0.0028 GeV², and the smearing radius σ_smear = 1.1047±0.0568 GeV. Compared with the original parameters, the optimized model improves the agreement with the 16 known charmonium states from the Particle Data Group. The root-mean-square deviation is reduced from 26.26 MeV to 24.62 MeV. Based on the optimized parameters, the mass spectrum of charmonium is predicted up to high excitation levels, including S, P, D, and F waves. The predicted energy level spacings and fine structures are consistent with the expected behavior of the confining potential, indicating that the optimized parameters retain the main physical features of the conventional potential model. Notably, the framework provides predictions for unobserved high-angular-momentum states such as the n³F_J and n¹F₃ series, which may serve as references for future experimental searches at facilities such as BESIII, Belle II and LHCb. The main innovation of this work is to combine neural network wavefunctions with statistical parameter inference, so that the solution of the Schrödinger equation, the fitting of potential parameters, and the quantification of parameter uncertainties are incorporated into an automated and extensible workflow. The proposed framework can be readily extended to other heavy quark systems such as bottomonium and B_c mesons, as well as to more complex spectroscopic problems involving coupled-channel effects.