This study uses the Marciniak and Kuczynski (M-K) method to present an analytical forming limit diagram (FLD) for sheet metals. The procedure for the analytical FLD prediction is described in detail and step-wise manner, and an algorithm is written using MATLAB. First, an appropriate algorithm is determined to establish the theoretical analyses, and various anisotropic yield functions, such as Hill’s 48, Barlat 89, and Hosford, are considered. The predicted FLDs are compared with experiments involving a typical AA6016-T4 aluminum alloy. Second, the Gurson model that considers damage growth is implemented when Hosford is the yield function, as Hosford criterion predicts the best comparable analytical FLD with experiments among the yield functions. Third, a parametric study is performed to investigate the effects of parameters on the FLD prediction. Results indicate that an extremely low value for the initial void volume fraction in the safe and groove zones has minimal effects on the FLD prediction. Lastly, the values of void volume fractions are calculated assuming no geometrical imperfections and the imperfection is because of higher void volume fraction in groove zone than that in safe zone.