In this paper, we consider the numerical approximation of high order Partial Differential Equations (PDEs) by means of NURBS-based Isogeometric Analysis (IGA) in the framework of the Galerkin method, for which global smooth basis functions with degree of continuity higher than C0 can be used. We derive a priori error estimates for high order elliptic PDEs under h-refinement, by extending existing results for second order PDEs approximated with IGA and specifically addressing the case of errors in lower order norms. We present some numerical results which both validate the proposed error estimates and highlight the accuracy of IGA. Then, we apply NURBS-based IGA to solve the fourth order stream function formulation of the Navier-Stokes equations for which we derive and numerically validate a priori error estimates under h-refinement. We solve the benchmark lid-driven cavity problem for Reynolds numbers up to 5000, by considering both the classical square cavity and the semi-circular cavity, which is exactly represented by NURBS.