We describe an efficient algorithm for the computation of separable isogenies between abelian varieties represented in the coordinate system given by algebraic theta functions. Let A be an abelian variety of dimension g defined over a field of odd characteristic. Our algorithm comprises two principal steps. First, given a theta null point for A and a subgroup K isotropic for the Weil pairing, we explain how to compute the theta null point corresponding to the quotient abelian variety A/K. Then, from the knowledge of a theta null point of A/K, we present an algorithm to obtain a rational expression for an isogeny from A to A/K. The algorithm that results from combining these two steps can be viewed as a higher-dimensional analog of the well-known algorithm of Vélu for computing isogenies between elliptic curves. In the case where K is isomorphic to (ℤ/ℓℤ)g for ℓ∈ℕ*, the overall time complexity of this algorithm is equivalent to O(log ℓ) additions in A and a constant number of ℓth root extractions in the base field of A. In order to improve the efficiency of our algorithms, we introduce a compressed representation that allows us to encode a point of level 4ℓ of a g-dimensional abelian variety using only g(g+1)/2⋅4g coordinates. We also give formulas for computing the Weil and commutator pairings given input points in theta coordinates.