An aerodynamic design method that treats multiple wings as well as their mutual interaction is proposed. This method provides section shapes of wings which realize specified surface pressure distributions. It can be applied to three dimensional design problems in potential, invisid and viscous, subsonic and transonic flow fields. The basic idea is an iterative residual reduction method using a new inverse problem solver which is the extension of Takanashi's integral equation method. An inverse problem is reformulated to be integral equations which express the relation of pressure differences to geometrical changes with consideration of interacting effects among wings. The inverse problem is solved numerically by introducing piecewise linear/constant function approximation. This method works well on several preliminary design problems. From both viewpoints of quality of solutions and efficiency, the method is found to be promising for complicated design problems.