The　soil　water　retention　curve　(SWRC)　is　essential　for　assessing　water　flow　and　solute　transport　in　unsaturated　media.　The　van　Genuchten　(VG)　model　is　widely　used　to　describe　the　SWRC;　however,　estimation　of　its　effective　hydraulic　parameters　is　often　prone　to　error,　especially　when　data　exist　for　only　a　limited　range　of　matric　potential.　We　developed　a　Metropolis-Hastings　algorithm　of　the　Markov　chain　Monte　Carlo　(MH-MCMC)　approach　using　R　to　estimate　VG　parameters,　which　produces　a　numerical　estimate　of　the　joint　posterior　distribution　of　model　parameters,　including　fully-quantified　uncertainties.　When　VG　model　parameters　were　obtained　using　complete　range　of　soil　water　content　(SWC)　data　(i.e.,　from　saturation　to　oven　dryness),　the　MH-MCMC　approach　returned　similar　accuracy　as　the　widely　used　non-linear　curve-fitting　program　RETC　(RETention　Curve),　but　avoiding　non-convergence　issues.　When　VG　model　parameters　were　obtained　using　5　SWC　data　measured　at　matric　potential　of　around　−60,　−100,　−200,　−500,　and　−15,000　cm,　the　MH-MCMC　approach　was　more　robust　than　the　RETC　program.　The　performance　of　MH-MCMC　are　generally　good　(R2　＞　0.95)　for　all　8　soils,　whereas　the　RETC　underperformed　for　coarse-textured　soils.　The　MH-MCMC　approach　was　used　to　obtain　VG　model　parameters　for　all　1871　soils　in　the　National　Cooperative　Soil　Characterization　dataset　with　SWC　measured　at　matric　potentials　of　−60　cm,　−100　cm,　−330　cm,　and　−15,000　cm;　the　results　showed　that　the　simulated　SWC　by　MH-MCMC　model　were　highly　consistent　with　the　measured　SWC　at　corresponding　matric　potential.　Altogether,　our　new　MH-MCMC　approach　to　solving　the　VG　model　is　more　robust　to　limited　coverage　of　soil　matric　potential　when　compared　to　the　RETC　procedures,　making　it　an　effective　alternative　to　traditional　water　retention　solvers.　We　developed　an　MH-MCMC　code　in　R　for　solving　VG　model　parameters,　which　can　be　found　at　the　GitHub　repository.　©　2022　by　the　authors.　Licensee　MDPI,　Basel,　Switzerland.