Bourgoyne and Young Method (BYM) is one of the most widely used rate of penetration (ROP) prediction methods. Drilling data, in this method, must be taken from uniform-lithology sections (homogeneous formations) to get the maximum prediction accuracy using multiple linear regression analysis. In the presence of heterogeneity, the accuracy of BYM decreases substantially due to complexity in lithology. There are various studies on ROP prediction based on BYM. Since these studies considers uniform lithology only, none of them functions satisfactorily in heterogeneous formations. Implementing a different type of regression model as an alternative to multiple linear regression by modifying BYM to predict ROP in heterogeneous environment for well programming is the main purpose of this study. This is done by introducing several geomechanical parameters to the BYM to estimate a relationship between the rock mechanical properties and the heterogeneity. An offset drilling data is taken from a heterogeneous formation and converted into BYM's functions. The effect of feature selection, data size, and different outlier labeling approaches on ROP prediction accuracy in heterogeneous formations are investigated by applying several predictive techniques: multiple linear regression, support vector regression, and artificial neural networks. The prediction results are compared and examined from a statistical point of view. It is revealed that the BYM could not represent ROP effectively in a heterogeneous environment. As demonstrated for the first time, the ROP prediction accuracy in heterogeneous formations significantly increases when machine learning techniques are used together with additional features comprising heterogeneity.