We propose a subgradient-based method for finding the maximum feasible subsystem in a collection of closed sets with respect to a given closed set $C$ (MFS$_C$). In this method, we reformulate the MFS$_C$ problem as an $\ell_0$ optimization problem and construct a sequence of continuous optimization problems to approximate it. The objective of each approximation problem is the sum of the composition of a nonnegative nondecreasing continuously differentiable concave function with the squared distance function to a closed set. Although this objective function is nonsmooth in general, a subgradient can be obtained in terms of the projections onto the closed sets. Based on this observation, we adapt a subgradient projection method to solve these approximation problems. Unlike classical subgradient methods, the convergence (clustering to stationary points) of our subgradient method is guaranteed with a {\em nondiminishing stepsize} under mild assumptions. This allows us to further study the sequential convergence of the subgradient method under suitable Kurdyka-{\L}ojasiewicz assumptions. Finally, we illustrate our algorithm numerically for solving the MFS$_C$ problems on a collection of halfspaces and a collection of unions of halfspaces, respectively, with respect to the set of $s$-sparse vectors.