This is an extended version of the workshop paper [1], in which a new computational technique called quantitative discrete approximation has been introduced. The technique provides finite discrete approximation of continuous dynamical systems which is suitable especially for a significant class of biochemical dynamical systems. With decreasing granularity the approximation of behaviour between a discrete state and its successor converges to the behaviour of the original continuous system in the respective part of the phase space. This paper provides a detailed description of the method and algorithms solving the reachability problem in biochemical dynamical systems. The method is supplemented with heuristics for reducing the cardinality of the reachable state space. The algorithms are evaluated on six models (with numbers of variables ranging from 2 to 12).