In several applications, as e.g. for problems in Solid Mechanics, Fluid Dynamics, and Biology, Partial Differential Equations (PDEs) are defined on surfaces. The numerical approximation of the PDEs often involves an approximation of the computational domain, represented by the surface, which may possibly lead to significant errors. Since a wide range of surfaces of practical interest can be represented by B-splines or NURBS, we consider the numerical solution of the PDEs by means of Isogeometric Analysis, a numerical method based on the isoparametric concept for which the same basis functions used to represent the computational domain are then used for approximating the unknown solution field of the PDEs. In this work, we solve linear, nonlinear, and time dependent PDEs involving the second order Laplace-Beltrami or high-order operators on surfaces by means of NURBS-based Isogeometric Analysis in the Galerkin framework. We highlight and discuss the properties of the method in terms of accuracy, efficiency, and regularity of the NURBS basis functions by performing error and spectrum analyses.