A design optimization method based on three dimensional Euler equations is developed. A finite volume method is implemented to discretize the Euler equations. Newton's method is used to solve the discretized form of Euler equations. Newton's method requires the calculation of the Jacobian matrix which is the derivative of the residual vector with respect to the flux vector. Different upwind methods are used in the calculation of flux vectors. Numerical and analytical methods are utilized in the evaluation of Jacobian matrices. The efficiency and accuracy of the analytical and numerical Jacobian evaluations are compared. In order to improve the accuracy of numerical method, detailed error analyses are performed. The optimum finite difference perturbation magnitude that minimizes the error is searched. The computation time of numerical Jacobian evaluation is reduced by calculating the flux vectors with perturbed flow variables only in related cells. The performances of different sparse matrix solvers are also compared. The effects of errors in numerical Jacobians on the accuracy of sensitivities is analyzed. Results show that the finite-difference perturbation magnitude and computer precision are the most important parameters that affect the accuracy of numerical Jacobians. Approximately the same optimum perturbation magnitude enables the most accurate numerical flux Jacobian and sensitivity calculations.