  A Skyrme-type effective potential is determined to describe the interaction
between $\Lambda$ hyperons in nuclear medium. Experimental data of the binding
energies of the double-$\Lambda$ ($\Lambda\Lambda$) nuclei with mass numbers
$A=10$--$13$ are used to fit the parameters of the $\Lambda\Lambda$
interaction. As a result of the fitting, we obtain eight different sets of the
$\Lambda\Lambda$ interaction parameters, which reproduces the input data within
5\% deviation from the experimental data on average. The eight $\Lambda\Lambda$
interactions are plugged in the calculation of the heavier $\Lambda\Lambda$
nuclei and the neutron star equation of state to explore the issue of hyperon
puzzle. We found that the $\Lambda\Lambda$ interaction, specifically, p-wave
interaction makes the equation of state stiff enough that the maximum mass of
the neutron star can be as large as, or above $2\;M_\odot$.
