We propose an algorithm to compute upper and lower bounds for the star discrepancy of an arbitrary sequence of points in the s-dimensional unit cube. The method is based on a particular partition of the unit cube into subintervals and on a specialized procedure for orthogonal range counting. The cardinality of the partition depends on the dimension and on an accuracy parameter that has to be specified. We have implemented this method and here we present results of some computational experiments obtained with this implementation.