A systematic diagrammatic expansion for Gutzwiller-wave functions (DE-GWF) proposed very recently is used for the description of superconducting (SC) ground state in the two-dimensional square-lattice $t$-$J$ model with the hopping electron amplitudes $t$ (and $t$) between nearest (and next-nearest) neighbors. On the example of the SC state analysis we provide a detailed comparison of the method results with other approaches. Namely: (i) the truncated DE-GWF method reproduces the variational Monte Carlo (VMC) results; (ii) in the lowest (zeroth) order of the expansion the method can reproduce the analytical results of the standard Gutzwiller approximation (GA), as well as of the recently proposed grand-canonical Gutzwiller approximation (GCGA). We obtain important features of the SC state. First, the SC gap at the Fermi surface resembles a $d_{x^2-y^2}$-wave only for optimally- and overdoped system, being diminished in the antinodal regions for the underdoped case in a qualitative agreement with experiment. Corrections to the gap structure are shown to arise from the longer range of the real-space pairing. Second, the nodal Fermi velocity is almost constant as a function of doping and agrees semi-quantitatively with experimental results. Third, we compare the doping dependence of the gap magnitude with experimental data. Fourth, we analyze the $mathbf{k}$-space properties of the model: Fermi surface topology and effective dispersion. The DE-GWF method opens up new perspectives for studying strongly-correlated systems, as: (i) it works in the thermodynamic limit, (ii) is comparable in accuracy to VMC, and (iii) has numerical complexity comparable to GA (i.e., it provides the results much faster than the VMC approach).