We consider the numerical approximation of geometric Partial Differential Equations (PDEs) defined on surfaces in the 3D space. In particular, we focus on the geometric PDEs deriving from the minimization of an energy functional by L2L2-gradient flow. We analyze two energy functionals: the area, which leads to the mean curvature flow, a nonlinear second order PDE, and the Willmore energy, leading to the Willmore flow, a nonlinear fourth order PDE. We consider surfaces represented by single-patch tensor product NURBS and discretize the PDEs by means of NURBS-based Isogeometric Analysis in the framework of the Galerkin method. To approximate the high order geometric PDEs we use high order continuous NURBS basis functions. For the time discretization of the nonlinear geometric PDEs, we use Backward Differentiation Formulas (BDF) with extrapolation of the geometric quantities involved in the weak formulation of the problem; in this manner, we solve a linear problem at each time step. We report numerical results concerning the mean curvature and Willmore flows on different geometries of interest and we show the accuracy and efficiency of the proposed approximation scheme.