In this study, the sheet metal forming limits are presented both in strain and stress spaces. The computations are based on strain gradient theory of plasticity in conjunction with Marciniak and Kuczynski model. This approach introduces an internal length scale into conventional constitutive equations and takes into account the effects of deformation inhomogeneity and material softening. The nonlinear differential equation of the thickness variation is the second order ordinary type. The shooting method has been used with different solvers such as Backward Differentiation Formula (BDFs or Gear method), explicit Runge–Kutta (the Dormand–Prince pair), explicit Runge–Kutta (the Bogacki–Shampine pair), Adams–Bashforth–Moulton and Rosenbrock formula and Trapezoidal rule. The numerical results of these techniques are almost the same. However, major differences can be observed in the speed of calculations. It was found that the explicit Runge–Kutta (the Dormand–Prince pair) initial value problem solver has the best speed of calculation. The forming limit results obtained in this work have been compared with some published experimental data.