In this paper we study the lepton favor violating decay channels of the neutral Higgs bosons of the Minimal Supersymmetric Standard Model into a lepton and an antilepton of different flavor. We work in the context of the most general flavor mixing scenario in the slepton sector, in contrast to the minimal flavor violation assumption more frequently used. Our analytic computation is a one-loop diagrammatic one, but in contrast to the full one-loop computation which is usually referred to the physical slepton mass basis, we use here instead the Mass Insertion Approximation (MIA) which uses the electroweak interaction slepton basis and treats perturbatively the mass insertions changing slepton flavor. By performing an expansion in powers of the external momenta in the relevant form factors, we will be able to separate explicitly in the analytic results the leading non-decoupling (constant at asymptotically large sparticle masses) and the next to leading decoupling contributions (decreasing with the sparticle masses). Our final aim is to provide a set of simple analytic formulas for the form factors and the associated effective vertices, that we think may be very useful for future phenomenological studies of the lepton flavor violating Higgs boson decays, and for their comparison with data. The accuracy of the numerical results obtained with the MIA are also analyzed and discussed here in comparison with the full one-loop results. Our most optimistic numerical estimates for the three neutral Higgs boson decays channels into τ and µ leptons, searching for their maximum rates that are allowed by present constraints from τ → µγ data and beyond Standard Model Higgs boson searches at the LHC, are also included.