In certain inverse problems it is useful to be able to compute solutions which are, in some sense, as simple as possible. For example,k one may wish to compute solutions which are piecewise constant and with as few discontinuities as possible. Such solutions are suited to describe models, e.g., geological layers, where the coarse structure is more important than the fine structure. A natural generalization of piecewise constant functions is piecewise polynomial solutions. In this paper we present a new algorithm which is capable of computing solutions that are piecewise polynomials, without having to specify a priori the positions of the break points between the polynomial pieces.