This article describes how to use smooth solvers for simulation of a class of piecewise smooth systems of ordinary differential equations, called Filippov systems, with discontinuous vector fields. In these systems constrained motion along a discontinuity surface (so-called sliding) is possible and requires special treatment numerically. The introduced algorithms are based on an extension to Filippov's method to stabilise the sliding flow together with accurate detection of the entrance and exit of sliding regions. The methods are implemented in a general way in MATLAB and sufficient details are given to enable users to modify the code to run on arbitrary examples. Here, the method is used to compute the dynamics of three example systems, a dry-friction oscillator, a relay feedback system and a model of an oil well drill-string.