Stratification reduces the variance of sample estimates for population parameters by creating homogeneous strata. Often, surveyors stratify the population using the most convenient variables such as age, sex, region, etc. Such convenient methods often do not produce internally homogeneous strata, hence, the precision of the estimates of the variables of interest could be further improved. This research addresses this problem and proposes a method for single study variables that follow a particular distribution estimated from the available data. The stratification problem is formulated as a mathematical programming problem and solved by using a dynamic programming technique. Methods for several distributions such as uniform, weibull, gamma, normal, log-normal, exponential, right-triangular and pareto are presented. An R package (called StratifyR) for the open-source 'The R Project for Statistical Computing' is created and implemented to construct optimal stratification boundaries and calculate sample sizes under Neyman allocation. Results reveal that the proposed method is able to compute precise stratum boundaries that are quite efficient when compared to other established methods.