This paper extends for the first time B \(\ddot{\text {o}}\) hmer/Schaback’s nonlinear discretization theory (J. Comput. Appl. Math., 254(2013), pp. 204–219) to nonlinear systems of partial differential equations. A direct collocation method is used to solve the coupled nonlinear systems. There are no restrictions on the nonlinear operators; in particular, there is no ellipticity, compactness, or self-adjointness assumed. The trial spaces might include a large class of finite-dimensional approximation spaces such as radial basis functions spaces, spline functions spaces, local Lagrange functions spaces, spectral techniques, and finite element spaces. A general framework for proving the convergence of collocation methods for solving well-posed nonlinear operator equations is presented. The theoretical results cover error bounds and convergence rates. In order to show how the general theoretical framework can be set to work, we take the meshfree method for solving the Navier–Stokes equations as an illustration and derive specific convergence rates in Sobolev spaces. Several numerical examples are provided to deal with numerical efficiency.