In this paper we present a scheme for the numerical solution of one-dimensional stochastic differential equations (SDEs) whose drift belongs to a fractional Sobolev space of negative regularity (a subspace of Schwartz distributions). We obtain a rate of convergence in a suitable $L^1$-norm and we implement the scheme numerically. To the best of our knowledge this is the first paper to study (and implement) numerical solutions of SDEs whose drift lives in a space of distributions. As a byproduct we also obtain an estimate of the convergence rate for a numerical scheme applied to SDEs with drift in $L^p$-spaces with $p\in(1,\infty)$.