This thesis deals with high angle of attack behaviour of a generic delta wing model aircraft. A high angle of attack wind tunnel database has been generated for this aircraft and based upon the bifurcation analysis of the data and the results of extensive simulations, it has been shown in the thesis that the post stall behaviour of this aircraft is both unstable and unpredictable. Unpredictability of aircraft behaviour arises from the fact that the aircraft response is oscillatory and divergent; the aircraft state trajectories do not settle down to any stable limit set and very often exceed valid aerodynamic database limits. This unpredictability of behaviour raises a major difficulty in the design of a procedure to recover the aircraft to normal flight regime in case the aircraft stalls and departs accidentally. A new methodology has been presented in this thesis to recover such unstable aircraft. In this methodology, a nonlinear controller is first designed at high angles of attack. This controller is connected by the pilot after the departure of the aircraft and the controller drives the aircraft to a well-defined spin condition. Thus, the controller makes the post stall aircraft behaviour predictable. Then a set of automatic recovery inputs is designed to reduce aircraft rotations and to lower the angle of attack. The present aircraft model is unstable at low angle of attack flight conditions as well and therefore to stabilize the aircraft to a low angle of attack level flight, another controller is designed. The high angle of attack controller is disconnected and the low angle of attack controller is connected automatically during the recovery process. The entire methodology is tested using extensive non-linear six degree-of-freedom simulations and the efficacy of the technique is established. The nonlinear controller that stabilizes the aircraft to a spin condition is designed using feedback linearization. The stability of a closed loop system obtained using feedback linearization is determined by the stability of the zero dynamics of the open loop plant. It has been shown in literature that the eigenvalues of the linearized zero dynamics are the same as the transmission zeros of the linearized plant at the equilibrium point. It is also well known that the location of transmission zeros of a linear system can be changed by the choice of outputs. In this thesis it is shown that if it is possible to reassign the outputs, then the feedback linearization based design for a linear system becomes very similar to a controller design for eigenvalue assignment. This thesis presents a new two-step procedure to obtain a locally stable and optimally robust closed loop system using feedback linearization. In the first step of this procedure optimal locations of the transmission zeros are found and in the second step, optimal outputs are constructed to place the system transmission zeros at these locations. The same outputs can then be used to construct nonlinear feedback for the nonlinear system and the resultant closed loop system is guaranteed to be locally robustly stable. The high angle of attack controller is designed using this procedure and its performance is presented in the thesis. The stabilized spin equilibrium point of the closed loop system is also shown to have a large domain of attraction. Having designed a locally robust stabilizing controller, the thesis addresses the problem of the evaluation of robustness of the stability of the equilibrium point in a nonlinear framework. The thesis presents a general method to construct bounds on the additive perturbations of the system vector field over a large region in the domain of attraction of a stable equilibrium point using Lyapunov functions. If the system perturbations lie within these bounds, the system is guaranteed to be stable. The thesis first proposes a method to numerically construct a Lyapunov function over a large region in the domain of attraction. In this method a sequence of Lyapunov functions are constructed such that each function in the sequence gives a larger estimate of the domain of attraction than the previous one. The seminal idea for this method is obtained from the existing literature and this idea is considerably generalized. Using this method, it is possible to numerically obtain a Lyapunov function value at each point in the domain of attraction, but the Lyapunov function does not have an analytical form. Hence, it is proposed to represent this function using neural networks. The thesis then discusses a new method to construct perturbation bounds. It is shown that the perturbation bounds obtained over a large region in the domain of attraction using a single Lyapunov function is too conservative. Using the concept of sequence of Lyapunov functions, the thesis proposes three methods to obtain the least conservative bounds for an initial local Lyapunov function. These general ideas are then applied to the aircraft example and the bounds on the perturbation of the aerodynamic database are presented.