**Abstract** : In this paper, we present a numerical algorithm for computing ISS Lyapunov func-tions for systems which are input-to-state stable (ISS) on compact subsets of the state space. The algorithm relies on a linear programming problem and computes a continuous piecewise affine ISS Lyapunov function on a simplicial grid covering the given compact set excluding a small neighbour-hood of the origin. We show that for every system admitting a twice continuously differentiable ISS Lyapunov function there exists a grid such that our algorithm terminates succesfully. 1. Introduction. The concept of input-to-state stability (ISS) was first intro-duced by Sontag [18] in the late 1980s and has soon turned out to be one of the most influential concepts for characterizing robust stability of nonlinear systems. Several results on ISS can be found in [18, 19, 20]. In [22], different equivalent formulations of ISS are given. In particular, it is shown that the ISS property is equivalent to the existence of an ISS Lyapunov function. The ISS notion is also useful in stability analysis of large scale systems. Stability of large interconnected networks of systems can be analyzed by means of ISS small gain theorems if all the subsystems are ISS [6, 7, 8, 14]. Motivated by these results, in this paper we focus on the computation of ISS Lyapunov functions for small systems, as the knowledge of ISS Lyapunov func-tions immediately leads to the knowledge of ISS gains which may be used in a small gain based stability analysis. Based on [21, Lemma 2.10–2.14], it was shown in [4] that ISS Lyapunov functions in implication form may be calculated using a Zubov approach. An alternative Zubov type approach was developed in [15]. In these two papers the problem of comput-ing ISS Lyapunov functions was transformed into the problem of computing robust Lyapunov functions for suitably designed auxiliary systems. This robust Lyapunov function can be characterized by the Zubov equation, a Hamilton-Jacobi-Bellman type partial differential equation, which can be solved numerically [3]. However, this approach only yields a numerical approximation of a Lyapunov function but not a true Lyapunov function. For discrete time systems, following the same auxiliary system approach, true Lyapunov functions can be computed by a set oriented approach [10]. This numerical approach, however, does not carry over to the continuous time setting. Moreover, the detour via the auxiliary system introduces conservatism, since the re-sulting Lyapunov function and ISS gains strongly depend on the way the auxiliary system is constructed.