  Numerical solving differential equations with fractional derivatives requires
elimination of the singularity which is inherent in the standard definition of
fractional derivatives. The method of integration by parts to eliminate this
singularity is well known. It allows to solve some equations but increases the
order of the equation and sometimes leads to wrong numerical results or
instability.
  We suggest another approach: the elimination of singularity by substitution.
It does not increase the order of equation and its numerical implementation
provides the opportunity to define fractional derivative as the limit of
discretization.
  We present a sufficient condition for the substitution-generated difference
approximation to be well-conditioned.
  We demonstrate how some equations can be solved using this method with full
confidence that the solution is accurate with at least second order of
approximation.
