Ordinary Differential Equations (ODEs) with discontinuous right-hand sides arise in a number of applications in elementary mechanics in relation to dry friction problems; in control theory in relation to switching controls and also in solving optimal control problems; in problems involving elastic-plastic behaviour and other variational inequalities. The crucial difference between these problems and problems with discontinuous coefficients is that the discontinuity is in the "state" variables, rather than occuring in the "time" variable. To emphasise this distinction we consider autonomous problems 1 = () where is discontinuous.The conventional theory regarding existence still works satisfactorily as long as the trajectory () is at a discontinuity for only finitely many . However, there are many examples where the trajectory is forced to stay at the discontinuity. In such a case conventional theory and conventional numerical methods in general break down; in particular, the existence of local solutions cannot be guaranteed (e.g. 1 = — sgn() + α if 0 < |α| < 1). The appropriate formulation of such a problem is to replace the ODE "1 = ()" by a differential inclusion 1 ∈ () where is an appropriate set-valued function with nonempty compact convex values and is upper semi-continuous (in terms of the "⊆" ordering). This formulation is due to Filippov(1960). With this formulation solutions can readily be shown to exist, although they are not necessarily unique.Numerical methods to compute solutions in the sense of Filippov have already been developed and convergence results proven by Taubert, Elliott and Niepage & Wendt. However, these methods are no better than 1st order accurate, and 1st order accuracy is established only with additional assumptions regarding the number of "switching points". Further, it is only guaranteed that a sub-sequence of the numerical solutions converge, unless solutions are unique. This limits the useful of these methods in solving problems where the desired solution is determined by other considerations, such as in two-point boundary value problems (TPBVPs).The approach taken here is to assume that the discontinuous ODE has some particular structure, and to exploit that to obtain high order accuracy and other information. The structure assumed is that there is a collection of disjoint open sets with dense union in , 1,..., such that for 1 ∈ implies that () = () with smooth. Further it is assumed that the can be described by indicator functions ℎ: = { | ℎ() < ℎ(), ≠ }. For each there is the active set () = { | ∈ ̄}. The Filippov formulation then becomes1 ∈ co{()| ∈ ()}.As determining the from is not a well posed problem we suppose that (()) is a piecewise constant function. On the intervals on which (()) is constant and given, we can select elements of the above differential inclusion that keep the active set constant and moreover gives a smooth ODE to solve. At a switching point the set of possible new values for (()) are determined by means of a Linear Complementarity Problem. Given a new active set selected from the list of possibilities, we can continue by solving a new smooth ODE to maintain the new active set, until another switching point occurs. Switching points can be both determined by the ODE itself, or if they are not unique, by the user of the algorithm.A complete algorithm for solving discontinuous ODEs by this approach has been developed and convergence results proven. These convergence results are of two types: Firstly, given only that the algorithm has produced a sequence of numerical solutions with step size approaching zero, then there is a converging subsequence and every converging subsequence converges to a solution. Secondly, if a particular solution is known which has only finitely many switching points and some non-degeneracy conditions are satisfied, then by suitable selections (in the algorithm) of new active sets a sequence of numerical solutions can be generated that converge to the given solution with the same order of accuracy as that of the smooth ODE solver used.This algorithm has been implemented and numerical results obtained. These results confirm the expectation of high order accuracy of this method, and that it is competitive with previously published methods.A number of extensions are proposed for this algorithm. The first is a "de-composition" type of approach that is particularly useful for friction problems with multiple friction surfaces. Convergence results are given for this extension. This extension has been implemented and numerical results obtained for a problem involving three friction surfaces.Some further extensions are also proposed to avoid some of the non-degeneracy conditions; these are in effect "higher order" modifications to the main algorithm. Also, some difficulties in applying the main algorithm to Hamiltonian inclusions (which are closely related to optimal control theory) are noted.