We investigate the practicability of an optimization algorithm based on the Broyden-Fletcher-Goldfarb-Shanno (BFGS) method for 3D shape design problems, using approximate sensitivities of objective functions, where the contribution of the partial derivatives of the flow state with respect to the control variables is neglected. Therefore, it is worthwhile to investigate how an optimization method based on the Hessian behaves in this context. Indeed, the Hessian should be far from its real value if the gradient approximation is wrong. The optimization methodology is characterized by an unstructured CAD-free framework for shape and mesh deformations, an automatic differentiation of programs for the computation of the gradient of the cost function, and an unstructured flow solver. The redesign of transonic and supersonic wings has been considered and the performance of the BFGS method has been analyzed in comparison with a steepest descent method. Taking into account that a line search is too expensive to be carried out in such problems, a step size proportional to the gradient modulus has been employed for updating the control variables. Numerical results show that the BFGS method does not suffer from the approximation used in the evaluation of sensitivities, and leads to an effective improvement of the efficiency of the optimization methodology. These results can be then considered an aposteriorijustification for incomplete sensitivities.