To be able to harness more power from the wind, wind turbine blades are getting longer. As they get longer, they get more flexible. This creates issues that have until recently not been of concern. Long flexible wind turbine blades can lose their stability to flow induced instabilities such as coupled-mode flutter. This type of flutter occurs when increasing wind speed causes a coupling of a bending and a torsional mode, which create limit cycle oscillations that can lead to blade failure. To be able to make the design of larger blades possible, it is important to be able to predict the critical flutter and post critical flutter behaviors of wind turbine blades. Most numerical research concerning coupled-mode wind turbine is focused on predicting the critical flutter point, and less focused on the post critical behavior. This is because of the mathematical complexities associated with the coupled, nonlinear wind turbine blade systems. Here, a numerical model is presented that predicts the critical flutter velocity and post critical flutter behavior for 3D airfoils with third order structural nonlinearities. The numerical model can account for the attached flow and separated flow region by using the ONERA dynamic stall model. By retaining higher-order structural nonlinearities, lateral and torsional displacements can be predicted, which makes it possible to use this model in the future to control wind turbine blade flutter. Furthermore, by using a dynamic stall model to simulate the flow, the solver is able to predict accurate limit cycle oscillations when the effective angle of attack is larger than the stall angle. The coupled, nonlinear equations of motion are two coupled nonlinear PDEs and are determined using Hamilton’s principle. In order to solve the equations of motion, they are discretized using the Galerkin technique into a set of ODEs. The motion of the airfoil is used as an input to ONERA. The airfoil is sectioned with the lateral position and angle of attack known as well as the velocity and acceleration of the section at an instance of time. This information is used by ONERA to calculate lift and moment coefficients for each section which are then used to calculate the total lift and moment forces of the airfoil. Then, a Fortran code solves the system by using Houbolt’s finite difference method. A theoretical NACA 0012 airfoil has been designed to define the parameters used by the equations of motion. Third bending and first torsional coupling occurs after the critical flutter point and dynamic lift and moment coefficients were observed. Dynamic stall was also observed at wind velocities farther away from the bifurcation point. Bifurcation diagrams, time histories, and phase planes have been created that represent the flutter behavior.