This paper describes how to use smooth solvers for simulation of a class of piecewise smooth dynamical systems, called Filippov systems, with discontinuous vector fields. In these systems constrained motion along a discontinuity surface (so-called sliding) is possible and require special treatment numerically. The introduced algorithms are based on an extension to Filippov's method to stabilize 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 arbitraray 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 a oil well drill-string.