diff --git a/Makefile.in b/Makefile.in
index d87350671a..cf0572bc1a 100644
--- a/Makefile.in
+++ b/Makefile.in
@@ -203,7 +203,8 @@ HEADER_DIRS :=                                                              \
         acb_mat         acb_poly        acb_calc        acb_hypgeom         \
         arb_fmpz_poly   arb_fpwrap                                          \
         acb_dft         acb_elliptic    acb_modular     acb_dirichlet       \
-        acb_theta       dirichlet       bernoulli       hypgeom             \
+        acb_theta       acb_ode                                             \
+        dirichlet       bernoulli       hypgeom                             \
                                                                             \
         gr              gr_generic      gr_vec          gr_mat              \
         gr_poly         gr_mpoly        gr_special                          \
diff --git a/dev/check_examples.sh b/dev/check_examples.sh
index 527ac10bdd..e639169d23 100755
--- a/dev/check_examples.sh
+++ b/dev/check_examples.sh
@@ -169,6 +169,9 @@ then
 elif test "$1" = "hilbert_matrix_ca";
 then
     echo "hilbert_matrix_ca....SKIPPED"
+elif test "$1" = "ode";
+then
+    echo "ode....SKIPPED"
 elif test "$1" = "huge_expr";
 then
     echo "huge_expr....SKIPPED"
diff --git a/doc/source/acb_poly.rst b/doc/source/acb_poly.rst
index 66148a38a0..f1ac950e15 100644
--- a/doc/source/acb_poly.rst
+++ b/doc/source/acb_poly.rst
@@ -1159,3 +1159,23 @@ Root-finding
     If this function returns zero, then the signs of the imaginary parts
     are not known for certain, based on the accuracy of the inputs
     and the working precision *prec*.
+
+
+Vector functions
+-------------------------------------------------------------------------------
+
+.. function:: acb_poly_struct * _acb_poly_vec_init(slong n)
+
+.. function:: void _acb_poly_vec_clear(acb_poly_struct *vec, slong n)
+
+.. function:: void _acb_poly_vec_set(acb_poly_struct *dest, const acb_poly_struct *src, slong n)
+
+.. function:: void _acb_poly_vec_set_block(acb_poly_struct *dest, const acb_poly_struct *src, slong n, slong base, slong len)
+
+.. function:: void _acb_poly_vec_fit_length(acb_poly_struct *vec, slong n, slong len)
+
+.. function:: void _acb_poly_vec_set_length(acb_poly_struct *vec, slong n, slong len)
+
+.. function:: void _acb_poly_vec_normalise(acb_poly_struct *vec, slong n)
+
+.. function:: int _acb_poly_vec_overlaps(acb_poly_struct *vec1, acb_poly_struct *vec2, slong n)
diff --git a/doc/source/gr.rst b/doc/source/gr.rst
index 79992c8a4d..8ae1845b5d 100644
--- a/doc/source/gr.rst
+++ b/doc/source/gr.rst
@@ -175,8 +175,8 @@ Flags can be OR'ed and checked only at the top level of a computation
 to avoid complex control flow::
 
     status = GR_SUCCESS;
-    gr |= gr_add(res, a, b, ctx);
-    gr |= gr_pow_ui(res, res, 2, ctx);
+    status |= gr_add(res, a, b, ctx);
+    status |= gr_pow_ui(res, res, 2, ctx);
     ...
 
 If we do not care about recovering from *undefined*/*unknown* results,
diff --git a/examples/ode.c b/examples/ode.c
new file mode 100644
index 0000000000..bfd613b51c
--- /dev/null
+++ b/examples/ode.c
@@ -0,0 +1,302 @@
+#include "acb_types.h"
+#include "acb_poly.h"
+#include "acb_mat.h"
+#include "acb_ode.h"
+#include "gr.h"
+#include "gr_poly.h"
+
+
+void
+_acb_ode_sum_swap_mat_ordinary(
+        acb_mat_t mat,
+        const acb_ode_sum_context_struct * ctx, slong s)
+{
+    for (slong j = 0; j < ctx->nsols; j++)
+    {
+        FLINT_ASSERT(ctx->sol[j].nlogs <= 1);
+        if (ctx->sol[j].nlogs < 1)
+            continue;
+        for (slong i = 0; i < ctx->nder; i++)
+            acb_swap(acb_mat_entry(mat, i, j),
+                     (ctx->sol[j].sums + s)->coeffs + i);
+    }
+}
+
+/* TODO
+ * - series
+ * - singular
+ */
+
+void
+ordinary(void)
+{
+    acb_ode_sum_context_t ctx;
+
+    /* (x^2 + 1)*Dx^3 + 7*x
+     * = (x^2 + 1)*Tx^3 + (-3*x^2 - 3)*Tx^2 + (2*x^2 + 2)*Tx + 7*x^4 */
+    acb_ode_sum_context_init(ctx, 4, 2, 3, 3);
+    acb_poly_set_coeff_si(ctx->dop + 3, 2, 1);
+    acb_poly_set_coeff_si(ctx->dop + 3, 0, 1);
+    acb_poly_set_coeff_si(ctx->dop + 2, 2, -3);
+    acb_poly_set_coeff_si(ctx->dop + 2, 0, -3);
+    acb_poly_set_coeff_si(ctx->dop + 1, 2, 2);
+    acb_poly_set_coeff_si(ctx->dop + 1, 0, 2);
+    acb_poly_set_coeff_si(ctx->dop + 0, 4, 7);
+
+    /* ctx->flags |= acb_ode_WANT_SERIES; */
+
+    acb_ode_sum_ordinary(ctx);
+    acb_ode_sum_canonical_basis(ctx);
+
+    acb_set_d_d(ctx->pts, 0.25, 0.25);
+    acb_set_si(ctx->pts + 1, 0);
+    mag_set_d(&((ctx->pts + 1)->real.rad), 1e-8);
+
+    ctx->prec = 64;
+    ctx->sums_prec = 64;
+
+    acb_ode_sum_divconquer(ctx, 20);
+
+    acb_mat_t mat;
+    acb_mat_init(mat, ctx->nder, ctx->nsols);
+
+    for (slong i = 0; i < ctx->npts; i++)
+    {
+        _acb_ode_sum_swap_mat_ordinary(mat, ctx, i);
+        acb_mat_printd(mat, 8);
+        /* flint_printf("%{acb_mat}\n\n", mat); */
+    }
+
+    acb_ode_sum_context_clear(ctx);
+    acb_mat_clear(mat);
+}
+
+
+void
+series(void)
+{
+    acb_ode_sum_context_t ctx;
+
+    acb_ode_sum_context_init(ctx, 2, 0, 1, 1);
+    acb_poly_set_coeff_si(ctx->dop + 1, 0, 1);
+    acb_poly_set_coeff_si(ctx->dop + 0, 1, -1);
+
+    acb_ode_sum_ordinary(ctx);
+    acb_ode_sum_canonical_basis(ctx);
+    ctx->prec = 64;
+    ctx->flags |= acb_ode_WANT_SERIES;
+
+    slong len = 5;
+    acb_ode_sum_divconquer(ctx, len);
+    /* Clear the high part reserved for the residual. (In this special case, the
+     * high part is zero because block_length = 1.) */
+    acb_poly_truncate(ctx->sol[0].series, len);
+
+    flint_printf("%{acb_poly}\n", ctx->sol[0].series);
+
+    acb_ode_sum_context_clear(ctx);
+}
+
+
+void
+bessel_j0(void)
+{
+    acb_ode_sum_context_t ctx;
+
+    slong dop_order = 2;
+    acb_ode_sum_context_init(ctx, dop_order + 1, 1, dop_order, dop_order);
+    acb_poly_set_coeff_si(ctx->dop + 2, 0, 1);
+    acb_poly_set_coeff_si(ctx->dop + 0, 2, 1);
+
+    ctx->sing_shifts[0].n = 0;
+    ctx->sing_shifts[0].mult = 2;
+
+    acb_ode_sum_canonical_basis(ctx);
+
+    acb_set_d(ctx->pts, 0.25);
+
+    ctx->prec = 64;
+    ctx->sums_prec = 64;
+
+    acb_ode_sum_divconquer(ctx, 10);
+
+    for (slong m = 0; m < ctx->nsols; m++)
+    {
+        acb_ode_sol_struct * sol = ctx->sol + m;
+        /* series in x */
+        flint_printf("f%wd =", m);
+        for (slong k = 0; k < sol->nlogs; k++)
+        {
+            flint_printf(" + (%{acb_poly})*log(%{acb} + x)^%wd/%wd",
+                         acb_ode_sol_sum_ptr(sol, 0, k),
+                         ctx->pts, k, k);
+        }
+        flint_printf("\n");
+    }
+
+    acb_ode_sum_context_clear(ctx);
+}
+
+
+void
+whittaker_m(void)  /* non-integer exponent, no logs */
+{
+    acb_t kappa, mu2, half;
+    acb_poly_t val;
+
+    acb_init(kappa);
+    acb_init(mu2);
+    acb_init(half);
+    acb_poly_init(val);
+
+    acb_ode_sum_context_t ctx;
+
+    acb_set_si(kappa, 2);
+    acb_set_si(mu2, 3);
+    acb_set_d(half, .5);
+
+    slong prec = 64;
+    slong dop_order = 2;
+    slong len = dop_order;
+
+    acb_ode_sum_context_init(ctx, dop_order + 1, 1, 1, len);
+
+    acb_poly_set_coeff_si(ctx->dop + 2, 0, 4);
+    acb_poly_set_coeff_si(ctx->dop + 1, 0, -4);
+    acb_poly_set_coeff_si(ctx->dop + 0, 2, -1);
+    acb_mul_si((ctx->dop + 0)->coeffs + 1, kappa, 4, prec);
+    acb_poly_set_coeff_si(ctx->dop + 0, 0, 1);
+    acb_addmul_si((ctx->dop + 0)->coeffs + 0, mu2, -4, prec);
+
+    acb_sqrt(ctx->expo, mu2, prec);
+    acb_add(ctx->expo, ctx->expo, half, prec);  /* other expo = mu2 - 1/2 */
+
+    ctx->sing_shifts[0].n = 0;
+    ctx->sing_shifts[0].mult = 1;
+
+    acb_ode_sum_canonical_basis(ctx);
+
+    acb_set_d(ctx->pts, 1.4242);
+
+    ctx->prec = prec;
+    ctx->sums_prec = prec;
+
+    acb_ode_sum_divconquer(ctx, prec);
+
+    acb_poly_struct * f = acb_ode_sol_sum_ptr(ctx->sol, 0, 0);
+
+    flint_printf("(%{acb} + x)^(%{acb}) * (%{acb_poly})\n",
+                 ctx->pts, ctx->expo, f);
+
+    _acb_ode_sol_value(val, ctx->expo, f, ctx->sol[0].nlogs, ctx->pts,
+                       ctx->nder, 1, ctx->prec);
+
+    flint_printf("M(%{acb} + x) = %{acb_poly} + O(x^%wd)\n", ctx->pts, val, len);
+
+    acb_ode_sum_context_clear(ctx);
+    acb_clear(kappa);
+    acb_clear(mu2);
+    acb_clear(half);
+    acb_poly_clear(val);
+}
+
+
+void
+fundamental_matrix(const char * dop_str,
+                   const acb_ode_exponents_struct * expos,
+                   double pt_d)
+{
+    gr_ctx_t CC, Pol, Dop;
+    gr_ptr dop;
+
+    slong prec = 30;
+
+    int status = GR_SUCCESS;
+
+    gr_ctx_init_complex_acb(CC, prec);
+    gr_ctx_init_gr_poly(Pol, CC);
+    gr_ctx_init_gr_poly(Dop, Pol);  /* should be Ore poly */
+
+    GR_TMP_INIT(dop, Dop);
+
+    status |= gr_ctx_set_gen_name(Pol, "z");
+    status |= gr_ctx_set_gen_name(Dop, "Tz");
+    status |= gr_set_str(dop, dop_str, Dop);
+    status |= gr_println(dop, Dop);
+    GR_MUST_SUCCEED(status);
+
+    slong dop_order = gr_poly_length(dop, Dop) - 1;
+
+    acb_mat_t mat;
+    acb_mat_init(mat, dop_order, dop_order);
+
+    acb_t pt;
+    acb_init(pt);
+    acb_set_d(pt, pt_d);
+
+    acb_ode_fundamental_matrix(mat, dop, Dop, expos, pt, 1, 0, 8, prec);
+
+    flint_printf("%{acb_mat}\n", mat);
+
+    acb_mat_clear(mat);
+    GR_TMP_CLEAR(dop, Dop);
+    gr_ctx_clear(Dop);
+    gr_ctx_clear(Pol);
+    gr_ctx_clear(CC);
+    acb_clear(pt);
+}
+
+
+void
+apery(void)
+{
+    acb_ode_shift_struct shift[1] = {{ .n = 0, .mult = 3 }};
+    acb_ode_group_struct grp[1] = {{ .nshifts = 1, .shifts = shift }};
+    acb_init(grp->expo);
+    acb_zero(grp->expo);
+    acb_ode_exponents_struct expos[1] = {{ .len = 1, .grps = grp }};
+
+    fundamental_matrix(
+            "(z^2 - 34*z + 1)*Tz^3 + (3*z^2 - 51*z)*Tz^2 + (3*z^2 - 27*z)*Tz + z^2 - 5*z",
+            expos,
+            0.015625);
+
+    acb_clear(grp->expo);
+}
+
+
+void
+multiple_shifts(void)
+{
+    /* const char * dop = "Tz^4 - 4*Tz^3 + 3*Tz^2 - z"; */
+    const char * dop = "Tz^6 - 6*Tz^5 + 12*Tz^4 - 10*Tz^3 + 3*Tz^2 + z^2";
+
+    acb_ode_shift_struct shift[3] = {
+        { .n = 0, .mult = 2 },
+        { .n = 1, .mult = 3 },
+        { .n = 3, .mult = 1 },
+    };
+    acb_ode_group_struct grp[1] = {{ .nshifts = 3, .shifts = shift }};
+    acb_init(grp->expo);
+    acb_zero(grp->expo);
+    acb_ode_exponents_struct expos[1] = {{ .len = 1, .grps = grp }};
+
+    fundamental_matrix(dop, expos, 2.);
+
+    acb_clear(grp->expo);
+}
+
+
+int
+main(void)
+{
+    /* ordinary(); */
+    /* series(); */
+    /* bessel_j0(); */
+    /* whittaker_m(); */
+
+    /* apery(); */
+    multiple_shifts();
+
+    flint_cleanup_master();
+}
diff --git a/src/acb_ode.h b/src/acb_ode.h
new file mode 100644
index 0000000000..6e4cded106
--- /dev/null
+++ b/src/acb_ode.h
@@ -0,0 +1,252 @@
+#ifndef acb_ode_H
+#define acb_ode_H
+
+/* temporary (for sage) */
+#include "gmp.h"
+typedef mp_limb_signed_t slong;
+/* */
+
+#include "acb_types.h"
+#include "gr_types.h"
+
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+        /****************************************************************************/
+
+        typedef struct
+        {
+                slong n;
+                slong mult;
+        }
+        acb_ode_shift_struct;
+
+
+        typedef struct
+        {
+                acb_t expo;
+                slong nshifts;
+                acb_ode_shift_struct * shifts;
+        }
+        acb_ode_group_struct;
+
+
+        typedef struct
+        {
+                slong len;
+                acb_ode_group_struct * grps;
+        }
+        acb_ode_exponents_struct;
+
+
+        typedef enum
+        {
+                acb_ode_BASIS_ECHELON = 0,
+                acb_ode_BASIS_CASCADE  /* XXX: same as Frobenius? */
+        }
+        acb_ode_basis_t;
+
+
+        slong acb_ode_group_length(const acb_ode_group_struct * grp);
+
+        /****************************************************************************/
+
+        typedef struct
+        {
+                acb_mat_t extini;
+
+                slong nlogs;
+                /* vector of polynomials in x, entries = coefficients of log(x)^k/k! of
+                 * a chunk of solution and/or rhs */
+                acb_poly_struct * series;
+                /* vector of polynomials holding the jets of coefficients wrt log(ξ)^k/k!
+                 * of the partial sums:
+                 * sums + j*alloc_logs + k is the jet of the coeff of log^k/k! at the
+                 * evaluation point of index j */
+                /* XXX switch to acb_struct*? we are not really using the polynomial
+                 * structure */
+                acb_poly_struct * sums;
+
+                slong alloc_logs;  /* XXX support automatic growth? */
+                slong alloc_pts;
+        }
+        acb_ode_sol_struct;
+
+        typedef acb_ode_sol_struct acb_ode_sol_t[1];
+
+        void acb_ode_sol_init(acb_ode_sol_struct * sol, slong nshifts,
+                              slong nlogs, slong npts);
+
+        void acb_ode_sol_clear(acb_ode_sol_struct * sol);
+
+        void acb_ode_sol_reset(acb_ode_sol_struct * sol);
+
+        void acb_ode_sol_fit_length(acb_ode_sol_struct * sol, slong len,
+                                    slong nder);
+
+        void acb_ode_sol_unit_ini(acb_ode_sol_struct * sol, slong i0,
+                                  const acb_ode_shift_struct * shifts);
+
+        acb_poly_struct * acb_ode_sol_sum_ptr(
+                const acb_ode_sol_struct * sol, slong j, slong k);
+
+        void
+                _acb_ode_sol_add_term(
+                        acb_ode_sol_struct * sol, slong base, slong n,
+                        acb_poly_struct * rhs, slong rhs_nlogs,
+                        /* XXX define a separate struct with the next three? */
+                        slong mult, slong ini_i, const acb_poly_struct * ind_n,
+                        acb_srcptr pows, const char * shifted, slong npts,
+                        const fmpz * binom_n,
+                        slong nder, slong prec, slong sums_prec);
+
+        void _acb_ode_sol_value(acb_poly_struct * val, acb_srcptr expo,
+                                const acb_poly_struct * sums, slong nlogs,
+                                acb_srcptr pt, slong nder,
+                                slong nshifts, slong prec);
+
+        /****************************************************************************/
+
+        /* XXX see if this can be reused by other summation algorithms (probably not as
+         * is), rename */
+
+#define acb_ode_WANT_SERIES 1
+
+        typedef struct
+        {
+                /* operator */
+                acb_poly_struct * dop;  /* XXX FJ: move outside the struct? */
+                slong dop_len;
+                slong dop_clen;
+
+                /* solution group */
+                acb_t expo;
+                acb_ode_shift_struct * sing_shifts;
+                acb_poly_t ind;
+
+                /* solutions */
+                acb_ode_sol_struct * sol;
+                slong nsols;
+
+                /* evaluation points */
+                acb_ptr pts;   /* evaluation points x[j] for 0 < j < npts */
+                acb_ptr pows;  /* x[j]^{n-i} for i < nder (i-major, cyclic on i);
+                                  only the block i=0 is really used for nonzero points */
+                char * shifted_sums;
+                slong npts;
+                slong nder;
+
+                /* working precisions */
+                slong prec;
+                slong sums_prec;
+                slong bounds_prec;
+
+                /* TODO: error bounds etc. */
+
+                /* ancillary data */
+                fmpz * binom_n;  /* binom(n, j) for j < nder */
+
+                /* options */
+                slong block_size;
+                ulong flags;
+
+                int have_precomputed;
+        }
+        acb_ode_sum_context_struct;
+
+        typedef acb_ode_sum_context_struct acb_ode_sum_context_t[1];
+
+        void acb_ode_sum_context_init(
+                acb_ode_sum_context_struct * ctx,
+                slong dop_len, slong npts, slong nsols, slong nder);
+        void acb_ode_sum_context_clear(acb_ode_sum_context_struct * ctx);
+
+        void acb_ode_sum_ordinary(acb_ode_sum_context_struct * ctx);
+        void acb_ode_sum_canonical_basis(acb_ode_sum_context_struct * ctx);
+        void acb_ode_sum_highest(acb_ode_sum_context_struct * ctx);
+        void acb_ode_sum_group(acb_ode_sum_context_struct * ctx, const acb_ode_group_struct * grp);
+
+
+        void _acb_ode_sum_precompute(acb_ode_sum_context_struct * ctx);
+        void _acb_ode_sum_reset(acb_ode_sum_context_struct * ctx);
+
+        void acb_ode_sum_divconquer(
+                acb_ode_sum_context_struct * ctx, slong nterms);
+
+        void acb_ode_sum_value(acb_poly_struct * val, slong nfrobshifts,
+                               const acb_ode_sum_context_struct * ctx,
+                               slong i, slong j, slong prec);
+
+        /****************************************************************************/
+
+        /* XXX void * data or fix args */
+        typedef void (* acb_ode_sum_worker_t)(
+                acb_ode_sum_context_struct * ctx, slong nterms);
+
+        void _acb_ode_fundamental_matrix(
+                acb_mat_struct * mat,
+                const acb_poly_struct * dop, slong dop_len,
+                const acb_ode_group_struct * groups, slong ngroups,
+                acb_srcptr pts, slong npts,
+                acb_ode_basis_t basis,
+                acb_ode_sum_worker_t sum_worker,
+                slong nterms, slong prec);
+
+        int acb_ode_fundamental_matrix(
+                acb_mat_struct * mat,
+                gr_srcptr dop, gr_ctx_t dop_ctx,
+                const acb_ode_exponents_struct * expos,
+                acb_srcptr pts, slong npts,
+                acb_ode_basis_t basis,
+                slong nterms, slong prec);
+
+        /****************************************************************************/
+
+        void _acb_ode_apply_diffop_basecase_weights(
+                acb_ptr weights,
+                const acb_poly_struct * dop, slong dop_len,
+                acb_srcptr expo, slong offset,
+                slong flen, slong nlogs, slong start, slong len,
+                slong prec);
+
+        void _acb_ode_apply_diffop_basecase_precomp(
+                acb_poly_struct * g, slong goff,
+                acb_srcptr weights, slong weights_nlogs,
+                const acb_poly_struct * f, slong foff, slong flen,
+                slong nlogs,
+                slong start, slong len,
+                slong prec);
+
+        void acb_ode_apply_diffop_basecase(
+                acb_poly_struct * g,
+                const acb_poly_struct * dop, slong dop_len,
+                acb_srcptr expo, slong offset,
+                const acb_poly_struct * f,
+                slong nlogs,
+                slong start, slong len,
+                slong prec);
+
+        void _acb_ode_apply_diffop_polmul(
+                acb_poly_struct * g, slong goff,
+                const acb_poly_struct * dop, slong dop_len,
+                acb_srcptr expo, slong offset,
+                const acb_poly_struct * f, slong foff, slong flen,
+                slong nlogs,
+                slong start, slong len,
+                slong prec);
+
+        void acb_ode_apply_diffop_polmul(
+                acb_poly_struct * g,
+                const acb_poly_struct * dop, slong dop_len,
+                acb_srcptr expo, slong offset,
+                const acb_poly_struct * f,
+                slong nlogs,
+                slong start, slong len,
+                slong prec);
+
+#ifdef __cplusplus
+}
+#endif
+
+#endif
diff --git a/src/acb_ode/apply_diffop_basecase.c b/src/acb_ode/apply_diffop_basecase.c
new file mode 100644
index 0000000000..499824f836
--- /dev/null
+++ b/src/acb_ode/apply_diffop_basecase.c
@@ -0,0 +1,164 @@
+#include "acb_types.h"
+#include "acb.h"
+#include "acb_poly.h"
+#include "fmpz_mat.h"
+#include "gr.h"
+#include "gr_mat.h"
+
+
+/* Notation:
+ *
+ * We are considering an operator
+ *
+ *   \sum_{i,j} a[i,j] x^j (x d/dx)^i
+ *
+ * applied to a series
+ *
+ *   f = x^{expo + offset} \sum_{p,k} f[k,p] x^p log(x)^k/k!,
+ *
+ * resulting in
+ *
+ *   g = x^{expo + offset} \sum{m1,k1} g[k1,m1] x^m1 log(x)^k1/k1!.
+ *
+ * We want to compute the terms of index m1 = start + p1 of g for 0 <= p1 < len.
+ */
+
+
+/* TODO specialized version for dop with fmpz coefficients (and expo == 0) */
+
+
+static void
+acb_addmul_binom(acb_ptr c, acb_srcptr b, fmpz_mat_t binom, slong i, slong t,
+                 slong prec)
+{
+    if (t == 0)
+        acb_add(c, c, b, prec);
+    else if (t == 1)
+        acb_addmul_si(c, b, i, prec);
+    else
+        acb_addmul_fmpz(c, b, fmpz_mat_entry(binom, i, t), prec);
+}
+
+
+void
+_acb_ode_apply_diffop_basecase_weights(
+        acb_ptr weights,  /* len * nlogs * flen */
+        const acb_poly_struct * dop, slong dop_len,
+        acb_srcptr expo, slong offset,
+        slong flen, slong nlogs, slong start, slong len,
+        slong prec)
+{
+
+    fmpz_mat_t binom;
+    if (nlogs >= 2)  /* XXX: worth caching? */
+    {
+        gr_ctx_t fmpz_ctx;
+        gr_ctx_init_fmpz(fmpz_ctx);
+        fmpz_mat_init(binom, dop_len, dop_len);
+        GR_MUST_SUCCEED(gr_mat_pascal((gr_mat_struct *) binom, -1, fmpz_ctx));
+        gr_ctx_clear(fmpz_ctx);
+    }
+
+    acb_t expo1;
+    acb_init(expo1);
+
+    if (expo != NULL && acb_is_zero(expo))
+        expo = NULL;
+
+    for (slong p1 = 0; p1 < len; p1++)  /* n1 = offset + start + p1 */
+    {
+        for (slong t = 0; t < nlogs; t++)  /* t = k - k1 */
+        {
+            for (slong p = 0; p < flen; p++)
+            {
+                slong j = start + p1 - p;
+                slong n = offset + p;
+
+                if (j < 0)
+                    break;
+
+                if (expo != NULL)
+                    acb_add_si(expo1, expo, n, prec);
+
+                acb_ptr c = weights + p1 * nlogs * flen + t * flen + p;
+                acb_zero(c);
+                for (slong i = dop_len - 1; i >= t; i--)  /* Horner */
+                {
+                    if (expo == NULL)
+                        acb_mul_si(c, c, n, prec);
+                    else
+                        acb_mul(c, c, expo1, prec);
+                    if (j >= (dop + i)->length)
+                        continue;
+                    acb_ptr aij = (dop + i)->coeffs + j;
+                    acb_addmul_binom(c, aij, binom, i, t, prec);
+                }
+            }
+        }
+    }
+
+    if (nlogs >= 2)
+        fmpz_mat_clear(binom);
+}
+
+
+void
+_acb_ode_apply_diffop_basecase_precomp(
+        acb_poly_struct * g, slong goff,
+        acb_srcptr weights, slong weights_nlogs,
+        const acb_poly_struct * f, slong foff, slong flen,
+        slong nlogs,
+        slong start, slong len,
+        slong prec)
+{
+    for (slong p1 = 0; p1 < len; p1++)  /* n1 = offset + start + p1 */
+    {
+        for (slong k1 = 0; k1 < nlogs; k1++)
+        {
+            acb_ptr dest = (g + k1)->coeffs + goff + p1;
+
+            /* XXX integrate the loop on k in acb_dot? */
+            for (slong k = k1; k < nlogs; k++)
+            {
+                acb_srcptr src = (f + k)->coeffs + foff;
+                acb_srcptr cofac = weights + p1 * weights_nlogs * flen + (k - k1) * flen;
+                slong fklen = FLINT_MIN(flen, (f + k)->length - foff);
+                /* loop on p */
+                acb_dot(dest, dest, 0, cofac, 1, src, 1, fklen, prec);
+            }
+        }
+    }
+}
+
+
+void
+acb_ode_apply_diffop_basecase(
+        acb_poly_struct * g,
+        const acb_poly_struct * dop, slong dop_len,
+        acb_srcptr expo, slong offset,
+        const acb_poly_struct * f,
+        slong nlogs,
+        slong start, slong len,
+        slong prec)
+{
+    slong flen = 0;
+    for (int k = 0; k < nlogs; k++)
+        flen = FLINT_MAX(flen, (f + k)->length);
+    /* XXX could we take flen = FLINT_MIN(flen, len) or similar? */
+
+    acb_ptr weights = _acb_vec_init(len * nlogs * flen);
+
+    _acb_ode_apply_diffop_basecase_weights(
+            weights, dop, dop_len, expo, offset,
+            flen, nlogs, start, len, prec);
+
+    _acb_poly_vec_fit_length(g, nlogs, len);
+    _acb_ode_apply_diffop_basecase_precomp(
+            g, 0,
+            weights, nlogs, f, 0, flen,
+            nlogs, start, len, prec);
+    _acb_poly_vec_set_length(g, nlogs, len);
+    _acb_poly_vec_normalise(g, nlogs);
+
+    _acb_vec_clear(weights, len * nlogs * flen);
+}
diff --git a/src/acb_ode/apply_diffop_polmul.c b/src/acb_ode/apply_diffop_polmul.c
new file mode 100644
index 0000000000..a2d005925c
--- /dev/null
+++ b/src/acb_ode/apply_diffop_polmul.c
@@ -0,0 +1,104 @@
+#include "acb_types.h"
+#include "acb_poly.h"
+
+
+static void
+xD_logcoeff_inplace(acb_poly_struct * f, acb_srcptr expo, slong offset, slong k,
+                    slong nlogs, slong prec)
+{
+    acb_poly_t tmp;
+    acb_poly_init(tmp);
+
+    /* XXX handle rational expo more efficiently (-> gr?) */
+    if (expo != NULL)
+        acb_poly_scalar_mul(tmp, f + k, expo, prec);
+    for (slong j = 0; j < (f + k)->length; j++)
+        acb_mul_ui((f + k)->coeffs + j, (f + k)->coeffs + j, offset + j, prec);
+    acb_poly_add(f + k, f + k, tmp, prec);
+    if (k + 1 < nlogs)
+        acb_poly_add(f + k, f + k, f + k + 1, prec);
+
+    acb_poly_clear(tmp);
+}
+
+
+/* Computes the coefficients of x^start to x^{start+len-1} (inclusive) of the
+ * log-polynomial g such that dop(x^{expo+offset}·f) = x^{expo+offset}·g.
+ *
+ * expo may be NULL (treated as zero)
+ *
+ * The underscore version operates on the chunk of coefficients of length
+ * flen starting at offset foff in each of the components of the input vector,
+ * and _adds_ each component of g at offset goff in the corresponding component
+ * of the output vector.
+ */
+
+/* XXX Take acb_struct**s instead of acb_poly_struct*s? More generally, how can
+ * we make the interface more natural / consistent with FLINT? */
+
+void
+_acb_ode_apply_diffop_polmul(
+        acb_poly_struct * g, slong goff,
+        const acb_poly_struct * dop, slong dop_len,
+        acb_srcptr expo, slong offset,
+        const acb_poly_struct * f, slong foff, slong flen,
+        slong nlogs,
+        slong start, slong len,
+        slong prec)
+{
+    acb_poly_t tmp;
+    acb_poly_init(tmp);
+    acb_poly_fit_length(tmp, start + len);
+
+    /* XXX leave it to the caller to do this, and destroy f? */
+    acb_poly_struct * curder = _acb_poly_vec_init(nlogs);
+    _acb_poly_vec_set_block(curder, f, nlogs, foff, flen);
+
+    for (slong i = 0; i < dop_len; i++)
+    {
+        for (slong k = 0; k < nlogs; k++)
+        {
+            /* should be a mulmid in the typical case (and should accept a
+             * pretransformed first operand) */
+            acb_poly_mullow(tmp, dop + i, curder + k, start + len, prec);
+            acb_poly_shift_right(tmp, tmp, start);
+
+            FLINT_ASSERT (tmp->length <= len);
+            _acb_poly_add((g + k)->coeffs + goff,
+                          (g + k)->coeffs + goff, len,
+                          tmp->coeffs, tmp->length,
+                          prec);
+
+            /* curder[k] ← (x·d/dx(prev curder))[k] */
+            xD_logcoeff_inplace(curder, expo, offset, k, nlogs, prec);
+        }
+    }
+
+    _acb_poly_vec_clear(curder, nlogs);
+    acb_poly_clear(tmp);
+}
+
+
+/* XXX support aliasing */
+void
+acb_ode_apply_diffop_polmul(
+        acb_poly_struct * g,
+        const acb_poly_struct * dop, slong dop_len,
+        acb_srcptr expo, slong offset,
+        const acb_poly_struct * f,
+        slong nlogs,
+        slong start, slong len,
+        slong prec)
+{
+    _acb_poly_vec_fit_length(g, nlogs, len);
+    _acb_ode_apply_diffop_polmul(
+            g, 0,
+            dop, dop_len,
+            expo, offset,
+            f, 0, start + len,
+            nlogs,
+            start, len,
+            prec);
+    _acb_poly_vec_set_length(g, nlogs, len);
+    _acb_poly_vec_normalise(g, nlogs);
+}
diff --git a/src/acb_ode/exponents.c b/src/acb_ode/exponents.c
new file mode 100644
index 0000000000..466358f01d
--- /dev/null
+++ b/src/acb_ode/exponents.c
@@ -0,0 +1,32 @@
+#include "acb_ode.h"
+
+
+slong
+acb_ode_group_length(const acb_ode_group_struct * grp)
+{
+    slong len = 0;
+
+    for (slong s = 0; s < grp->nshifts; s++)
+        len += grp->shifts[s].mult;
+
+    return len;
+}
+
+
+void
+acb_ode_exponents_init(acb_ode_exponents_struct * expos)
+{
+    expos->len = 0;
+    expos->grps = NULL;
+}
+
+
+void
+acb_ode_exponents_free(acb_ode_exponents_struct * expos)
+{
+    flint_free(expos->grps);
+}
+
+
+
+
diff --git a/src/acb_ode/fundamental_matrix.c b/src/acb_ode/fundamental_matrix.c
new file mode 100644
index 0000000000..b8bc235809
--- /dev/null
+++ b/src/acb_ode/fundamental_matrix.c
@@ -0,0 +1,210 @@
+#include "acb_types.h"
+#include "acb_poly.h"
+#include "acb_mat.h"
+#include "acb_ode.h"
+
+#include "gr.h"
+#include "gr_poly.h"
+
+
+static void
+init_context(acb_ode_sum_context_t ctx,
+             const acb_poly_struct * dop, slong dop_len,
+             const acb_ode_group_struct * grp,
+             acb_srcptr pts, slong npts, slong nder,
+             slong prec)
+{
+    acb_ode_sum_context_init(ctx, dop_len, npts, grp->nshifts, nder);
+    _acb_poly_vec_set(ctx->dop, dop, dop_len);
+    acb_ode_sum_group(ctx, grp);
+    _acb_vec_set(ctx->pts, pts, npts);
+    ctx->prec = prec;        /* XXX */
+    ctx->sums_prec = prec;
+}
+
+
+static slong
+col(const acb_ode_group_struct * grp, slong s, slong k)
+{
+    slong j = 0;
+    for (slong s1 = 0; s1 <= s; s1++)
+        j += grp->shifts[s1].mult;
+    return j - 1 - k;
+}
+
+
+static void
+fill_column(acb_mat_struct * mat,
+            const acb_ode_group_struct * grp, slong s, slong k,
+            const acb_poly_struct * val)
+{
+    slong j = col(grp, s, k);
+    slong len = FLINT_MIN(acb_mat_nrows(mat), val->length);
+    for (slong i = 0; i < len; i++)
+    {
+        acb_ptr a = acb_mat_entry(mat, i, j);
+        acb_swap(a, val->coeffs + i);
+    }
+}
+
+
+static void
+fix_column_echelon(acb_mat_struct * mat,
+                   const acb_ode_group_struct * grp, slong s, slong k,
+                   const acb_mat_t extini, slong prec)
+{
+    slong mult = grp->shifts[s].mult;
+    slong delta = mult - 1 - k;
+    slong j = col(grp, s, k);
+
+    for (slong s1 = s + 1, j1 = col(grp, s1, 0); s1 < grp->nshifts; s1++)
+    {
+        slong mult1 =  grp->shifts[s1].mult;
+        for (slong k1 = FLINT_MAX(0, mult1 - delta); k1 < mult1; k1++, j1++)
+        {
+            acb_ptr cc = acb_mat_entry(extini, s1, k1 + delta);
+            /* flint_printf("fix: s=%wd k=%wd s1=%wd k1=%wd j=%wd j1=%d cc=%{acb}\n", s, k, s1, k1, j, j1, cc); */
+            for (slong i = 0; i < acb_mat_nrows(mat); i++)
+            {
+                acb_ptr a  = acb_mat_entry(mat, i, j);
+                acb_ptr a1 = acb_mat_entry(mat, i, j1);
+                /* flint_printf("i=%wd %{acb} - %{acb}", i, a, a1); */
+                acb_submul(a, cc, a1, prec);
+                /* flint_printf(" = %{acb}\n", a); */
+            }
+        }
+    }
+
+}
+
+
+static void
+fill_group(acb_mat_struct * mat, const acb_ode_sum_context_t ctx,
+           const acb_ode_group_struct * grp, slong p,
+           acb_ode_basis_t basis, slong prec)
+{
+    for (slong s = grp->nshifts - 1; s >= 0; s--)
+    {
+        acb_poly_struct * val = _acb_poly_vec_init(ctx->sol[s].nlogs);
+
+        slong mult = grp->shifts[s].mult;
+        acb_ode_sum_value(val, mult, ctx, s, p, prec);
+        /* flint_printf("s=%wd mult=%wd val=%{acb_poly}\n\n", s, mult, val); */
+
+        for (slong k = 0; k < mult; k++)
+        {
+            slong delta = mult - 1 - k;
+            fill_column(mat, grp, s, k, val + delta);
+
+            /* flint_printf("s=%wd k=%wd after fill \n%{acb_mat}\n\n", s, k, mat); */
+
+            switch (basis)
+            {
+                case acb_ode_BASIS_CASCADE:
+                    break;
+                case acb_ode_BASIS_ECHELON:
+                    fix_column_echelon(mat, grp, s, k, ctx->sol[s].extini, prec);
+                    break;
+                default:
+                    FLINT_ASSERT(0);
+            }
+
+            /* flint_printf("s=%wd k=%wd after fix \n%{acb_mat}\n\n", s, k, mat); */
+        }
+
+        _acb_poly_vec_clear(val, mult);
+    }
+}
+
+
+void
+_acb_ode_fundamental_matrix(
+        acb_mat_struct * mat,
+        const acb_poly_struct * dop, slong dop_len,
+        const acb_ode_group_struct * groups, slong ngroups,
+        acb_srcptr pts, slong npts,
+        acb_ode_basis_t basis,
+        acb_ode_sum_worker_t sum_worker,
+        slong nterms, slong prec)
+{
+    slong nder = acb_mat_nrows(mat);
+
+    for (slong g = 0, j = 0; g < ngroups; g++)
+    {
+        slong glen = acb_ode_group_length(groups + g);
+
+        acb_ode_sum_context_t ctx;
+        init_context(ctx, dop, dop_len, groups + g, pts, npts, nder, prec);
+        /* XXX tmp */
+        /* ctx->flags |= acb_ode_WANT_SERIES; */
+
+        sum_worker(ctx, nterms);
+
+        for (slong p = 0; p < npts; p++)
+        {
+            acb_mat_t win;
+            acb_mat_window_init(win, mat + p, j, 0, j + glen, acb_mat_nrows(mat));
+
+            fill_group(win, ctx, groups + g, p, basis, prec);
+
+            acb_mat_window_clear(win);
+        }
+
+        acb_ode_sum_context_clear(ctx);
+
+        j += glen;
+    }
+}
+
+
+int
+acb_ode_fundamental_matrix(
+        acb_mat_struct * mat,
+        gr_srcptr dop, gr_ctx_t dop_ctx,
+        const acb_ode_exponents_struct * expos,
+        acb_srcptr pts, slong npts,                   /* XXX gr_vec? */
+        acb_ode_basis_t basis,
+        slong nterms, slong prec)
+{
+    gr_ctx_t CC, Pol, Dop;
+    gr_ptr _dop;
+
+    if (dop_ctx->which_ring != GR_CTX_GR_POLY)
+        return GR_UNABLE;
+
+    ulong which_base = POLYNOMIAL_ELEM_CTX(dop_ctx)->which_ring;
+    if (!(which_base == GR_CTX_GR_POLY || which_base == GR_CTX_FMPZ_POLY ||
+          which_base == GR_CTX_FMPQ_POLY))
+        return GR_UNABLE;
+
+    int status = GR_SUCCESS;
+
+    gr_ctx_init_complex_acb(CC, prec);
+    gr_ctx_init_gr_poly(Pol, CC);
+    gr_ctx_init_gr_poly(Dop, Pol);  /* XXX should be Ore poly */
+
+    status |= gr_ctx_set_gen_name(Dop, POLYNOMIAL_CTX(dop_ctx)->var);
+    status |= gr_ctx_set_gen_name(
+            Pol, POLYNOMIAL_CTX(POLYNOMIAL_ELEM_CTX(dop_ctx))->var);
+
+    GR_TMP_INIT(_dop, Dop);
+    status |= gr_set_other(_dop, dop, dop_ctx, Dop);
+
+    GR_MUST_SUCCEED(status);
+
+    gr_ptr dop_coeff = gr_poly_entry_ptr(_dop, 0, Dop);
+    slong dop_len = gr_poly_length(_dop, Dop);
+
+    _acb_ode_fundamental_matrix(mat, dop_coeff, dop_len,
+                                expos->grps, expos->len,
+                                pts, npts,
+                                basis, acb_ode_sum_divconquer,
+                                nterms, prec);
+
+    GR_TMP_CLEAR(_dop, Dop);
+    gr_ctx_clear(Dop);
+    gr_ctx_clear(Pol);
+    gr_ctx_clear(CC);
+
+    return status;
+}
diff --git a/src/acb_ode/sol.c b/src/acb_ode/sol.c
new file mode 100644
index 0000000000..640d8feaf8
--- /dev/null
+++ b/src/acb_ode/sol.c
@@ -0,0 +1,96 @@
+#include "acb_mat.h"
+#include "acb_poly.h"
+#include "acb_ode.h"
+
+
+void
+acb_ode_sol_init(acb_ode_sol_struct * sol, slong nshifts,
+                 slong nlogs, slong npts)
+{
+    acb_mat_init(sol->extini, nshifts, nlogs);
+
+    sol->nlogs = 0;
+    sol->series = _acb_poly_vec_init(nlogs);
+
+    sol->sums = _acb_poly_vec_init(npts * nlogs);
+
+    sol->alloc_logs = nlogs;
+    sol->alloc_pts = npts;
+}
+
+
+void
+acb_ode_sol_set_ini(acb_ode_sol_struct * sol, acb_srcptr ini,
+                    const acb_ode_shift_struct * shifts)
+{
+    slong nshifts = acb_mat_nrows(sol->extini);
+    slong i = 0;
+    for (slong j = 0; j < nshifts; j++)
+    {
+        slong mult = shifts ? shifts[j].mult : 1;
+        for (slong k = 0; k < mult; k++)
+            acb_set(acb_mat_entry(sol->extini, j, k), ini + i++);
+    }
+}
+
+
+void
+acb_ode_sol_unit_ini(acb_ode_sol_struct * sol, slong i0,
+                     const acb_ode_shift_struct * shifts)
+{
+    slong nshifts = acb_mat_nrows(sol->extini);
+    for (slong i = 0, j = 0; j < nshifts; j++)
+    {
+        slong mult = shifts ? shifts[j].mult : 1;
+        for (slong k = 0; k < mult; k++, i++)
+            acb_set_si(acb_mat_entry(sol->extini, j, k), i == i0 ? 1 : 0);
+    }
+}
+
+
+void
+acb_ode_sol_clear(acb_ode_sol_struct * sol)
+{
+    acb_mat_clear(sol->extini);
+    _acb_poly_vec_clear(sol->sums, sol->alloc_pts * sol->alloc_logs);
+    _acb_poly_vec_clear(sol->series, sol->alloc_logs);
+}
+
+
+void
+acb_ode_sol_fit_length(acb_ode_sol_struct * sol, slong len,
+                       slong nder)
+{
+    for (slong k = 0; k < sol->alloc_logs; k++)
+    {
+        acb_poly_struct * f = sol->series + k;
+        acb_poly_fit_length(f, len);
+        _acb_poly_set_length(f, f->alloc);  /* for printing */
+    }
+
+    for (slong i = 0; i < sol->alloc_pts * sol->alloc_logs; i++)
+    {
+        acb_poly_struct * s = sol->sums + i;
+        acb_poly_fit_length(s, nder);
+        _acb_poly_set_length(s, s->alloc);  /* for printing */
+    }
+}
+
+
+void
+acb_ode_sol_reset(acb_ode_sol_struct * sol)
+{
+    for (slong k = 0; k < sol->alloc_logs; k++)
+        acb_poly_zero(sol->series + k);
+
+    for (slong i = 0; i < sol->alloc_pts * sol->alloc_logs; i++)
+        acb_poly_zero (sol->sums + i);
+}
+
+
+acb_poly_struct *
+acb_ode_sol_sum_ptr(const acb_ode_sol_struct * sol,
+                    slong j, slong k)
+{
+    return sol->sums + j * sol->alloc_logs + k;
+}
diff --git a/src/acb_ode/sol_next_term.c b/src/acb_ode/sol_next_term.c
new file mode 100644
index 0000000000..3b11fc46da
--- /dev/null
+++ b/src/acb_ode/sol_next_term.c
@@ -0,0 +1,163 @@
+#include "acb_types.h"
+#include "acb.h"
+#include "acb_mat.h"
+#include "acb_ode.h"
+
+
+typedef enum
+{
+    LC_UNKNOWN,
+    LC_ACB,
+    LC_FMPZ
+}
+lc_status_t;
+
+
+/* XXX take off = n - base instead of (base, n)?
+ * XXX separate rhs_off from off? */
+
+
+static void
+next_coeff(acb_ode_sol_struct * sol, slong base, slong n,
+           acb_poly_struct * rhs, slong rhs_nlogs,
+           slong mult, slong ini_i, const acb_poly_struct * ind_n,
+           slong prec)
+{
+    fmpz_t lc_fmpz;
+    acb_t invlc;
+    slong k;
+
+    lc_status_t lc_status = LC_UNKNOWN;
+    fmpz_init(lc_fmpz);
+    acb_init(invlc);
+
+    while (rhs_nlogs > 0
+           && acb_is_zero((rhs + rhs_nlogs - 1)->coeffs + n - base))
+        rhs_nlogs--;
+
+    /* Compute the high-log part of the new term from the coefficient of x^{λ+n}
+     * on the rhs and the indicial polynomial */
+
+    acb_struct * new_term = _acb_vec_init(rhs_nlogs);
+
+    for (k = rhs_nlogs - 1; k >= 0; k--)
+    {
+        /* combin = rhs[k][0] + sum(ind_n[mult+1+j] * new_term[k+1+j], j≥0) */
+        acb_dot(new_term + k,
+                (rhs + k)->coeffs + n - base,
+                0,
+                ind_n->coeffs + mult + 1, 1,
+                new_term + k + 1, 1,
+                rhs_nlogs - k - 1,
+                prec);
+
+        if (acb_is_zero(new_term + k))
+            continue;
+
+        /* only compute invlc if needed */
+        if (lc_status == 0)
+        {
+            if (acb_is_exact(ind_n->coeffs + mult)
+                && acb_get_unique_fmpz(lc_fmpz, ind_n->coeffs + mult))
+            {
+                lc_status = LC_FMPZ;
+                acb_indeterminate(invlc);
+            }
+            else
+            {
+                lc_status = LC_ACB;
+                acb_inv(invlc, ind_n->coeffs + mult, prec);
+            }
+        }
+
+        if (lc_status == LC_FMPZ)
+            acb_div_fmpz(new_term + k, new_term + k, lc_fmpz, prec);
+        else
+            acb_mul(new_term + k, new_term + k, invlc, prec);
+
+        acb_neg(new_term + k, new_term + k);
+    }
+
+    /* Write the new term to sol, store new extini */
+
+    for (k = 0; k < mult; k++)  /* initial conditions */
+    {
+        acb_set((sol->series + k)->coeffs +  n - base,
+                acb_mat_entry(sol->extini, ini_i, k));
+    }
+    for (; k < rhs_nlogs + mult; k++)
+    {
+        acb_swap((sol->series + k)->coeffs + n - base, new_term + k - mult);
+
+        if (mult > 0)
+            acb_set(acb_mat_entry(sol->extini, ini_i, k),
+                    (sol->series + k)->coeffs + n - base);
+    }
+    for (; k < sol->nlogs; k++)
+    {
+        acb_zero((sol->series + k)->coeffs + n - base);
+    }
+
+    /* Update log-degree and used initial values */
+
+    sol->nlogs = FLINT_MAX(sol->nlogs, rhs_nlogs + mult);
+
+    _acb_vec_clear(new_term, rhs_nlogs);
+    acb_clear(invlc);
+    fmpz_clear(lc_fmpz);
+}
+
+
+static void
+next_sums(acb_ode_sol_struct * sol, slong base, slong n,
+          acb_srcptr pows, const char * shifted, slong npts,
+          const fmpz * binom_n,
+          slong nder, slong prec)
+{
+    acb_t tmp;
+    acb_init(tmp);
+
+    for (slong j = 0; j < npts; j++)
+    {
+        for (slong k = 0; k < sol->nlogs; k++)
+        {
+            acb_ptr c = (sol->series + k)->coeffs + n - base;
+            acb_ptr s = acb_ode_sol_sum_ptr(sol, j, k)->coeffs;
+
+            if (shifted[j])
+            {
+                slong pows_offset = (n % nder) * npts;
+                acb_srcptr ptpow = pows + pows_offset + j;
+                acb_mul(tmp, c, ptpow, prec);
+                for (slong i = 0; i < nder; i++)
+                    acb_addmul_fmpz(s + i, tmp, binom_n + i, prec);
+            }
+            else
+            {
+                for (slong i = 0; i < nder; i++)
+                {
+                    slong pows_offset = ((n + nder - i) % nder) * npts;
+                    acb_srcptr ptpow = pows + pows_offset + j;
+                    acb_mul(tmp, c, ptpow, prec);
+                    acb_addmul_fmpz(s + i, tmp, binom_n + i, prec);
+                }
+            }
+        }
+    }
+
+    acb_clear(tmp);
+}
+
+
+void
+_acb_ode_sol_add_term(
+        acb_ode_sol_struct * sol, slong base, slong n,
+        acb_poly_struct * rhs, slong rhs_nlogs,
+        slong mult, slong ini_i, const acb_poly_struct * ind_n,
+        acb_srcptr pows, const char * shifted, slong npts,
+        const fmpz * binom_n,
+        slong nder, slong prec, slong sums_prec)
+{
+    next_coeff(sol, base, n, rhs, rhs_nlogs, mult, ini_i, ind_n, prec);
+    next_sums(sol, base, n, pows, shifted, npts, binom_n, nder, prec);
+}
diff --git a/src/acb_ode/sol_values.c b/src/acb_ode/sol_values.c
new file mode 100644
index 0000000000..3c3b9ec8c5
--- /dev/null
+++ b/src/acb_ode/sol_values.c
@@ -0,0 +1,73 @@
+#include "acb_types.h"
+#include "acb_poly.h"
+#include "acb_ode.h"
+
+
+void
+_acb_ode_sol_value(acb_poly_struct * val, acb_srcptr expo,
+                   const acb_poly_struct * sums, slong nlogs,
+                   acb_srcptr pt, slong nder,
+                   slong nfrobshifts, slong prec)
+{
+    acb_poly_t expt, exptpow, log, logpow, term;
+
+    acb_poly_init(expt);
+    acb_poly_init(exptpow);
+    acb_poly_init(log);
+    acb_poly_init(logpow);
+    acb_poly_init(term);
+
+    acb_poly_set_coeff_si(expt, 1, 1);
+    acb_poly_set_coeff_acb(expt, 0, pt);
+
+    /* XXX does this work when pt contains 0 but expo is a nonneg int? */
+    acb_poly_pow_acb_series(exptpow, expt, expo, nder, prec);
+    acb_poly_log_series(log, expt, nder, prec);
+
+    acb_poly_one(logpow);
+
+    if (nlogs >= 1)
+    {
+        slong len = FLINT_MIN(nlogs, nfrobshifts);
+        _acb_poly_vec_set(val, sums, len);
+        _acb_poly_vec_zero(val + len, nlogs - len);
+    }
+
+    for (slong p = 1; p < nlogs; p++)
+    {
+        acb_poly_mullow(logpow, logpow, log, nder, prec);
+        _acb_vec_scalar_div_ui(logpow->coeffs, logpow->coeffs, logpow->length,
+                               p, prec);
+
+        for (slong d = 0; d < FLINT_MIN(nlogs - p, nfrobshifts); d++)
+        {
+            acb_poly_mullow(term, sums + d + p, logpow, nder, prec);
+            acb_poly_add(val + d, val + d, term, prec);
+        }
+    }
+
+    for (slong d = 0; d < nfrobshifts; d++)
+        acb_poly_mullow(val + d, val + d, exptpow, nder, prec);
+
+    acb_poly_clear(expt);
+    acb_poly_clear(exptpow);
+    acb_poly_clear(log);
+    acb_poly_clear(logpow);
+    acb_poly_clear(term);
+}
+
+
+/* XXX rename? move? */
+
+void
+acb_ode_sum_value(acb_poly_struct * val, slong nfrobshifts,
+                  const acb_ode_sum_context_struct * ctx,
+                  slong i, slong j, slong prec)  /* sol #i, pt #j */
+{
+    /* flint_printf("i=%wd j=%wd nlogs=%wd\n", i, j, ctx->sol[i].nlogs); */
+    /* flint_printf("%{acb_poly}\n", acb_ode_sol_sum_ptr(ctx->sol + i, j, 0)); */
+    _acb_ode_sol_value(val, ctx->expo,
+                       acb_ode_sol_sum_ptr(ctx->sol + i, j, 0),
+                       ctx->sol[i].nlogs, ctx->pts + j, ctx->nder,
+                       nfrobshifts, prec);
+}
diff --git a/src/acb_ode/sum_context.c b/src/acb_ode/sum_context.c
new file mode 100644
index 0000000000..5744a67466
--- /dev/null
+++ b/src/acb_ode/sum_context.c
@@ -0,0 +1,162 @@
+#include <string.h>
+
+#include "acb_types.h"
+#include "acb.h"
+#include "acb_mat.h"
+#include "acb_poly.h"
+#include "acb_ode.h"
+#include "fmpz_vec.h"
+
+
+typedef acb_ode_sum_context_struct * ctx_ptr;
+
+
+void
+acb_ode_sum_context_init(ctx_ptr ctx, slong dop_len, slong npts,
+                         slong nsols, slong nder)
+{
+    slong dop_order = dop_len - 1;
+
+    ctx->dop = _acb_poly_vec_init(dop_len);
+    ctx->dop_len = dop_len;
+
+    acb_init(ctx->expo);
+
+    acb_poly_init(ctx->ind);
+
+    ctx->sing_shifts = flint_calloc(2 * dop_order, sizeof(slong));
+    for (slong s = 0; s < dop_order; s++)
+        ctx->sing_shifts[s].n = -1;
+
+    ctx->sol = flint_malloc(nsols * sizeof(acb_ode_sol_struct));
+    /* using dop_order as a bound
+     * - for max possible log prec (will need updating to support inhomogeneous
+     *   equations),
+     * - for number of initial value positions */
+    for (slong m = 0; m < nsols; m++)
+        acb_ode_sol_init(ctx->sol + m, dop_order, dop_order, npts);
+    ctx->nsols = nsols;
+
+    ctx->pts = _acb_vec_init(npts);
+    ctx->pows = _acb_vec_init(npts * nder);
+    ctx->shifted_sums = flint_malloc(npts);
+    ctx->npts = npts;
+    ctx->nder = nder;
+
+    ctx->bounds_prec = MAG_BITS;  /* default value */
+
+    ctx->binom_n = _fmpz_vec_init(nder);
+
+    ctx->flags = 0;
+    ctx->have_precomputed = 0;
+}
+
+
+void
+acb_ode_sum_context_clear(ctx_ptr ctx)
+{
+    _acb_poly_vec_clear(ctx->dop, ctx->dop_len);
+
+    acb_clear(ctx->expo);
+    acb_poly_clear(ctx->ind);
+
+    flint_free(ctx->sing_shifts);
+
+    for (slong m = 0; m < ctx->nsols; m++)
+        acb_ode_sol_clear(ctx->sol + m);
+    flint_free(ctx->sol);
+
+    _acb_vec_clear(ctx->pts, ctx->npts);
+    _acb_vec_clear(ctx->pows, ctx->npts * ctx->nder);
+    flint_free(ctx->shifted_sums);
+
+    _fmpz_vec_clear(ctx->binom_n, ctx->nder);
+}
+
+
+void
+acb_ode_sum_ordinary(ctx_ptr ctx)
+{
+    for (slong n = 0; n < ctx->dop_len - 1; n++)
+    {
+        ctx->sing_shifts[n].n = n;
+        ctx->sing_shifts[n].mult = 1;
+    }
+}
+
+
+void
+acb_ode_sum_canonical_basis(ctx_ptr ctx)
+{
+    for (int m = 0; m < ctx->nsols; m++)
+        acb_ode_sol_unit_ini(ctx->sol + m, m, ctx->sing_shifts);
+}
+
+
+void
+acb_ode_sum_highest(ctx_ptr ctx)
+{
+    for (int m = 0; m < ctx->nsols; m++)
+    {
+        acb_mat_zero(ctx->sol[m].extini);
+        acb_one(acb_mat_entry(ctx->sol[m].extini, m,
+                              ctx->sing_shifts[m].mult - 1));
+    }
+}
+
+
+void
+acb_ode_sum_group(ctx_ptr ctx, const acb_ode_group_struct * grp)
+{
+    acb_set(ctx->expo, grp->expo);
+    memcpy(ctx->sing_shifts, grp->shifts,
+           grp->nshifts * sizeof(acb_ode_shift_struct));
+    acb_ode_sum_highest(ctx);
+}
+
+
+void
+_acb_ode_sum_precompute(ctx_ptr ctx)
+{
+    if (ctx->have_precomputed)
+        return;
+
+    slong dop_len = ctx->dop_len;
+
+    ctx->dop_clen = 0;
+    for (slong i = 0; i < dop_len; i++)
+        ctx->dop_clen = FLINT_MAX(ctx->dop_clen, ctx->dop[i].length);
+
+    ctx->block_size = FLINT_MAX(1, ctx->dop_clen - 1);
+
+    acb_poly_fit_length(ctx->ind, dop_len);
+    _acb_poly_set_length(ctx->ind, dop_len);
+    for (slong i = 0; i < dop_len; i++)
+        acb_poly_get_coeff_acb((ctx->ind)->coeffs + i, ctx->dop + i, 0);
+    _acb_poly_normalise(ctx->ind);
+
+    /* We can trade some full-width muls for muls by small integers by
+     * multiplying everyone by ξ^n), but this optimization is only
+     * valid when ξ does not contain 0. */
+    for (slong j = 0; j < ctx->npts; j++)  /* XXX move? */
+        ctx->shifted_sums[j] = !acb_contains_zero(ctx->pts + j);
+
+    /* TODO error bounds etc. */
+
+    ctx->have_precomputed = 1;
+}
+
+
+void
+_acb_ode_sum_reset(ctx_ptr ctx)
+{
+    for (slong m = 0; m < ctx->nsols; m++)
+    {
+        acb_ode_sol_reset(ctx->sol + m);
+        acb_ode_sol_fit_length(ctx->sol + m, 2 * ctx->block_size,
+                               ctx->nder);
+    }
+
+    for (slong i = 0; i < ctx->npts; i++)
+        acb_one(ctx->pows + i);
+}
diff --git a/src/acb_ode/sum_divconquer.c b/src/acb_ode/sum_divconquer.c
new file mode 100644
index 0000000000..0751accb2e
--- /dev/null
+++ b/src/acb_ode/sum_divconquer.c
@@ -0,0 +1,337 @@
+#include "acb_types.h"
+#include "acb.h"
+#include "acb_poly.h"
+#include "acb_ode.h"
+#include "fmpz.h"
+#include "fmpz_vec.h"
+
+
+static void
+acb_poly_shift_series(acb_poly_struct * g, const acb_poly_struct * f,
+                      acb_srcptr a, slong n, slong len, slong prec)
+{
+    acb_poly_fit_length(g, len);
+
+    if (len == 1 && acb_is_zero(a))  /* covers ordinary points */
+    {
+        acb_ptr c = g->coeffs;
+        acb_zero(c);
+
+        for (slong i = acb_poly_degree(f); i >= 0; i--)
+        {
+            acb_mul_si(c, c, n, prec);
+            acb_add(c, c, f->coeffs + i, prec);
+        }
+
+        _acb_poly_set_length(g, 1);
+        _acb_poly_normalise(g);
+    }
+    else
+    {
+        acb_t b;
+        acb_init(b);
+
+        acb_add_si(b, a, n, prec);
+        acb_poly_taylor_shift(g, f, b, prec);  /* wasteful */
+        acb_poly_truncate(g, len);
+
+        acb_clear(b);
+    }
+}
+
+
+static slong
+max_nlogs(acb_ode_sum_context_struct * ctx)
+{
+    slong nlogs = 0;
+    for (slong m = 0; m < ctx->nsols; m++)
+        nlogs = FLINT_MAX(nlogs, ctx->sol[m].nlogs);
+    return nlogs;
+}
+
+
+static slong
+max_nlogs_xn(acb_ode_sum_context_struct * ctx, slong off)
+{
+    for (slong nlogs = max_nlogs(ctx); nlogs > 0; nlogs--)
+    {
+        for (slong m = 0; m < ctx->nsols; m++)
+        {
+            if (nlogs <= ctx->sol[m].nlogs)
+            {
+                acb_poly_struct * p = ctx->sol[m].series + nlogs - 1;
+                if (!acb_is_zero(p->coeffs + off))
+                    return nlogs;
+            }
+        }
+    }
+    return 0;
+}
+
+
+/* XXX better decouple rhs from sol */
+
+
+static void
+apply_diffop(acb_ode_sum_context_struct * ctx,
+             slong base, slong low, slong mid, slong high)
+{
+    /* flint_printf("apply_diffop: base=%wd low=%wd mid=%wd high=%wd\n"); */
+
+    if (high - mid >= 7 * ctx->dop_len)  /* threshold from ore_algebra */
+    {
+        for (slong m = 0; m < ctx->nsols; m++)
+            _acb_ode_apply_diffop_polmul(
+                    ctx->sol[m].series, mid - base,
+                    ctx->dop, ctx->dop_len,
+                    ctx->expo, low,
+                    ctx->sol[m].series, low - base, mid - low,
+                    ctx->sol[m].nlogs,
+                    mid - low, high - mid,
+                    ctx->prec);
+    }
+    else
+    {
+        slong nlogs = max_nlogs(ctx);
+
+        slong weights_len = (high - mid) * nlogs * (mid - low);
+        acb_ptr weights = _acb_vec_init(weights_len);
+
+        _acb_ode_apply_diffop_basecase_weights(
+                weights, ctx->dop, ctx->dop_len, ctx->expo,
+                low, mid - low, nlogs, mid - low, high - mid,
+                ctx->prec);
+
+        for (slong m = 0; m < ctx->nsols; m++)
+            _acb_ode_apply_diffop_basecase_precomp(
+                    ctx->sol[m].series, mid - base,
+                    weights, nlogs, ctx->sol[m].series, low - base, mid - low,
+                    ctx->sol[m].nlogs, mid - low, high - mid, ctx->prec);
+
+        _acb_vec_clear(weights, weights_len);
+    }
+}
+
+
+static void
+next_binom_n(fmpz * binom_n, slong n, slong len)
+{
+    if (n == 0)
+    {
+        _fmpz_vec_zero(binom_n, len);
+        fmpz_one(binom_n);
+    }
+    else
+    {
+        for (slong i = len - 1; i > 0; i--)
+            fmpz_add(binom_n + i, binom_n + i, binom_n + i - 1);
+    }
+}
+
+
+void
+_acb_ode_sum_forward_1(acb_ode_sum_context_struct * ctx,
+                       slong base, slong n)
+{
+    slong ini_i, mult, len;
+    acb_poly_t ind_n;
+
+    acb_poly_init(ind_n);
+
+    FLINT_ASSERT(n >= base);
+
+    /* find the position and multiplicity of the current shift as an indicial
+     * root */
+
+    mult = 0;
+    for (slong i = 0; i < ctx->dop_len - 1; i++)
+    {
+        /* flint_printf("i=%wd shift=%wd mult=%wd\n", i, ctx->sing_shifts[i].n, ctx->sing_shifts[i].mult); */
+        if (ctx->sing_shifts[i].n == n)
+        {
+            ini_i = i;
+            mult = ctx->sing_shifts[i].mult;
+            break;
+        }
+    }
+
+    /* find the maximum log(x)-length of the corresponding coefficient of x on
+     * the rhs */
+
+    len = max_nlogs_xn(ctx, n - base);
+
+    /* evaluate the indicial polynomial */
+
+    acb_poly_shift_series(ind_n, ctx->ind, ctx->expo, n, len + mult, ctx->prec);
+
+    next_binom_n(ctx->binom_n, n, ctx->nder);
+
+    /* flint_printf("fwd1: n=%wd mult=%wd len=%wd ind_n=%{acb_poly}\n", n, mult, len, ind_n); */
+
+    /* advance solutions */
+
+    for (slong m = 0; m < ctx->nsols; m++)
+    {
+        acb_ode_sol_struct * sol = ctx->sol + m;
+        _acb_ode_sol_add_term(
+                sol, base, n,
+                sol->series, sol->nlogs, mult, ini_i, ind_n,
+                ctx->pows, ctx->shifted_sums, ctx->npts, ctx->binom_n, ctx->nder,
+                ctx->prec, ctx->sums_prec);
+    }
+
+    /* compute the next power of each evaluation point */
+    /* XXX wasteful near the end */
+
+    for (slong j = 0; j < ctx->npts; j++)
+    {
+        acb_mul(ctx->pows + ((n + 1) % ctx->nder) * ctx->npts + j,
+                ctx->pows + (n       % ctx->nder) * ctx->npts + j,
+                ctx->pts + j,
+                ctx->sums_prec);
+    }
+
+    acb_poly_clear(ind_n);
+}
+
+
+static void
+discard_block(acb_ode_sum_context_struct * ctx, slong size)
+{
+    for (slong m = 0; m < ctx->nsols; m++)
+    {
+        for (slong k = 0; k < ctx->sol[m].nlogs; k++)
+        {
+            acb_poly_struct * f = ctx->sol[m].series + k;
+            /* NOTE the sage version uses slightly different sizes near
+             * max_terms */
+            _acb_poly_shift_right(f->coeffs, f->coeffs, 2 * size, size);
+            _acb_vec_zero(f->coeffs + size, size);  /* XXX useful? */
+        }
+    }
+}
+
+
+void
+_acb_ode_sum_divconquer(acb_ode_sum_context_struct * ctx,
+                        slong base, slong low, slong high)
+{
+    slong mid;
+
+    slong bs = ctx->block_size;
+
+    FLINT_ASSERT(base <= low);
+    FLINT_ASSERT(low <= high);
+
+    if (high == low)
+        return;
+    else if (high == low + 1) {
+        _acb_ode_sum_forward_1(ctx, base, low);
+        return;
+    }
+
+    if (high <= low + bs)
+        mid = (high + low)/2;
+    else
+        mid = low + bs;
+
+    _acb_ode_sum_divconquer(ctx, base, low, mid);
+
+    /*
+       acb_poly_t tmp;
+       acb_poly_init(tmp);
+       flint_printf("divconquer: bs=%wd base=%wd low=%wd mid=%wd high=%wd\n", bs, base, low, mid, high);
+       for (slong m = 0; m < ctx->nsols; m++)
+       {
+       acb_poly_set(tmp, ctx->sol[m].series);
+       acb_poly_truncate(tmp, mid - base);
+       flint_printf("series[%wd]=%{acb_poly}\n", m, tmp);
+       }
+       acb_poly_clear(tmp);
+       */
+
+    slong resid_len = FLINT_MIN(high - mid, bs);
+
+    apply_diffop(ctx, base, low, mid, mid + resid_len);
+
+    /*
+       for (slong m = 0; m < ctx->nsols; m++)
+       flint_printf("series+res[%wd]=%{acb_poly}\n", m, ctx->sol[m].series);
+       for (slong j = 0; j < ctx->npts; j++)
+       for (slong m = 0; m < ctx->nsols; m++)
+       flint_printf("sums[%wd][%wd]=%{acb_poly}\n", m, j, ctx->sol[m].sums + j);
+       */
+
+    /* TODO convergence check (requires the residual!) */
+
+    /* TODO adjust working precisions */
+
+    if (mid - low >= bs)
+    {
+        if(ctx->flags & acb_ode_WANT_SERIES)
+        {
+            for (slong m = 0; m < ctx->nsols; m++)
+                acb_ode_sol_fit_length(ctx->sol + m, high + bs, 0);
+        }
+        else
+        {
+            discard_block(ctx, bs);
+            base += bs;
+        }
+    }
+
+    /* XXX in WANT_SERIES mode, set the correct lengths before terminating? */
+
+    /* XXX We sometimes do not compute the residual corresponding to the last
+     * block. This may need to change when adding support for error bounds. */
+    _acb_ode_sum_divconquer(ctx, base, mid, high);
+}
+
+
+void
+_acb_ode_sum_fix(acb_ode_sum_context_struct * ctx)
+{
+    acb_t inv, invpow;
+    acb_init(inv);
+    acb_init(invpow);
+
+    for (slong j = 0; j < ctx->npts; j++)
+    {
+        if (!ctx->shifted_sums[j])
+            continue;
+
+        acb_inv(inv, ctx->pts + j, ctx->sums_prec);
+
+        acb_one(invpow);
+        for (slong i = 1; i < ctx->nder; i++)
+        {
+            acb_mul(invpow, invpow, inv, ctx->sums_prec);
+            for (slong m = 0; m < ctx->nsols; m++)
+            {
+                for (slong k = 0; k < ctx->sol[m].nlogs; k++)
+                {
+                    acb_poly_struct * f = acb_ode_sol_sum_ptr(
+                            ctx->sol + m, j, k);
+                    acb_ptr c = f->coeffs + i;
+                    acb_mul(c, c, invpow, ctx->sums_prec);
+                }
+            }
+        }
+    }
+
+    acb_clear(invpow);
+    acb_clear(inv);
+}
+
+
+/* XXX temporary, to be generalized */
+void
+acb_ode_sum_divconquer(
+        acb_ode_sum_context_struct * ctx,
+        slong nterms)
+{
+    _acb_ode_sum_precompute(ctx);
+    _acb_ode_sum_reset(ctx);
+    _acb_ode_sum_divconquer(ctx, 0, 0, nterms);
+    _acb_ode_sum_fix(ctx);
+}
diff --git a/src/acb_ode/test/main.c b/src/acb_ode/test/main.c
new file mode 100644
index 0000000000..1f49478779
--- /dev/null
+++ b/src/acb_ode/test/main.c
@@ -0,0 +1,10 @@
+#include "t-apply_diffop.c"
+#include "t-sum_divconquer.c"
+
+test_struct tests[] =
+{
+    TEST_FUNCTION(acb_ode_apply_diffop),
+    TEST_FUNCTION(acb_ode_sum_divconquer)
+};
+
+TEST_MAIN(tests)
diff --git a/src/acb_ode/test/t-apply_diffop.c b/src/acb_ode/test/t-apply_diffop.c
new file mode 100644
index 0000000000..cad8281871
--- /dev/null
+++ b/src/acb_ode/test/t-apply_diffop.c
@@ -0,0 +1,99 @@
+#include "test_helpers.h"
+#include "acb_types.h"
+#include "acb_poly.h"
+#include "fmpq_types.h"
+#include "fmpq_poly.h"
+#include "acb_ode.h"
+
+TEST_FUNCTION_START(acb_ode_apply_diffop, state)
+{
+    fmpq_poly_t tmp;
+    fmpq_poly_init(tmp);
+
+    for (slong iter = 0; iter < 10000 * flint_test_multiplier(); iter++)
+    {
+        slong dop_len = n_randint(state, 3);
+        acb_poly_struct * dop = _acb_poly_vec_init(dop_len);
+        for (slong i = 0; i < dop_len; i++)
+        {
+            acb_poly_randtest(dop + i, state, 1 + n_randint(state, 10), 1 +
+                              n_randint(state, 200), 5);
+            /* fmpq_poly_randtest(tmp, state, 1 + n_randint(state, 5), 1 + n_randint(state, 5)); */
+            /* flint_printf("dop[%wd] = %{fmpq_poly}\n", i, tmp); */
+            /* acb_poly_set_fmpq_poly(dop + i, tmp, 100); */
+        }
+
+        acb_t expo;
+        acb_init(expo);
+        acb_randtest(expo, state, 1 + n_randint(state, 200), 5);
+        slong offset = n_randint(state, 1000);
+
+        slong logs = n_randint(state, 4);
+        acb_poly_struct * f = _acb_poly_vec_init(logs);
+        for (slong k = 0; k < logs; k++)
+        {
+            acb_poly_randtest(f + k, state, 1 + n_randint(state, 20), 1 +
+                              n_randint(state, 200), 5);
+            /* fmpq_poly_randtest(tmp, state, 1 + n_randint(state, 10), 1 + n_randint(state, 5)); */
+            /* flint_printf("f[%wd] = %{fmpq_poly}\n", k, tmp); */
+            /* acb_poly_set_fmpq_poly(f + k, tmp, 100); */
+        }
+
+        slong len1 = 1 + n_randint(state, 5);
+        slong len2 = n_randint(state, 5);
+        slong len = len1 + len2;
+        slong start = n_randint(state, len);
+
+        slong prec = 2 + n_randint(state, 100);
+
+        acb_poly_struct * g1 = _acb_poly_vec_init(logs);
+        acb_poly_struct * g2 = _acb_poly_vec_init(logs);
+        acb_poly_struct * g3 = _acb_poly_vec_init(logs);
+
+        acb_ode_apply_diffop_polmul(g1, dop, dop_len, expo, offset, f, logs, start, len, prec);
+
+        _acb_poly_vec_fit_length(g2, logs, len1 + len2);
+        /* TODO play with foff and offset too */
+        _acb_ode_apply_diffop_polmul(g2, 0, dop, dop_len, expo, offset, f, 0, start + len1, logs, start, len1, prec);
+        _acb_ode_apply_diffop_polmul(g2, len1, dop, dop_len, expo, offset, f, 0, start + len, logs, start + len1, len2, prec);
+        _acb_poly_vec_set_length(g2, logs, len);
+        _acb_poly_vec_normalise(g2, logs);
+
+        acb_ode_apply_diffop_basecase(g3, dop, dop_len, expo, offset, f, logs, start, len, prec);
+
+        if (!_acb_poly_vec_overlaps(g1, g2, logs)
+            || !_acb_poly_vec_overlaps(g1, g3, logs)
+            || !_acb_poly_vec_overlaps(g2, g3, logs))
+        {
+            flint_printf("FAIL\n\n");
+            for (slong i = 0; i < dop_len; i++)
+                flint_printf("dop[%wd] = %{acb_poly}\n", i, dop + i);
+            for (slong k = 0; k < logs; k++)
+                flint_printf("f[%wd] = %{acb_poly}\n", k, f + k);
+            flint_printf("expo = %{acb}\n", expo);
+            flint_printf("offset = %wd\n", offset);
+            flint_printf("start = %wd\n", start);
+            flint_printf("len1 = %wd\n", len1);
+            flint_printf("len2 = %wd\n", len2);
+            for (slong k = 0; k < logs; k++)
+            {
+                flint_printf("g1[%wd] = %{acb_poly}\n", k, g1 + k);
+                flint_printf("g2[%wd] = %{acb_poly}\n", k, g2 + k);
+                flint_printf("g3[%wd] = %{acb_poly}\n", k, g3 + k);
+            }
+            flint_abort();
+        }
+
+        /* multiplier à g/d par des trucs */
+        /* quelques séries connues ? */
+
+        // TODO aliasing tests
+
+        _acb_poly_vec_clear(dop, dop_len);
+        acb_clear(expo);
+        _acb_poly_vec_clear(g1, logs);
+        _acb_poly_vec_clear(g2, logs);
+    }
+
+    TEST_FUNCTION_END(state);
+}
diff --git a/src/acb_ode/test/t-sum_divconquer.c b/src/acb_ode/test/t-sum_divconquer.c
new file mode 100644
index 0000000000..79a4b9a09c
--- /dev/null
+++ b/src/acb_ode/test/t-sum_divconquer.c
@@ -0,0 +1,101 @@
+#include "test_helpers.h"
+#include "acb_types.h"
+#include "acb_poly.h"
+#include "fmpq_types.h"
+#include "fmpq_poly.h"
+#include "acb_ode.h"
+
+TEST_FUNCTION_START(acb_ode_sum_divconquer, state)
+{
+
+    for (slong iter = 0; iter < 10 * flint_test_multiplier(); iter++)
+    {
+        /* classical functions -- well for now just atan */
+
+        acb_ode_sum_context_t ctx;
+
+        slong dop_order = 2;
+        acb_ode_sum_context_init(ctx, dop_order + 1,
+                                 n_randint(state, 5),
+                                 dop_order,
+                                 1 + n_randint(state, 3));
+
+        acb_poly_set_coeff_si(ctx->dop + 2, 2, 1);
+        acb_poly_set_coeff_si(ctx->dop + 2, 0, 1);
+        acb_poly_set_coeff_si(ctx->dop + 1, 2, 1);
+        acb_poly_set_coeff_si(ctx->dop + 1, 0, -1);
+
+        acb_ode_sum_ordinary(ctx);
+        acb_ode_sum_canonical_basis(ctx);
+
+        ctx->prec = n_randint(state, 1000);
+        ctx->sums_prec = ctx->prec;
+
+        for (slong i = 0; i < ctx->npts; i++)
+        {
+            acb_ptr a = ctx->pts + i;
+            acb_randtest_precise(a, state, ctx->prec, 0);
+            acb_div_si(a, a, 4, ctx->prec);
+        }
+
+        ctx->flags = n_randint(state, 2);
+
+        slong nterms = ctx->prec;
+        acb_ode_sum_divconquer(ctx, nterms);
+
+        acb_poly_struct * ref = _acb_poly_vec_init(dop_order);
+
+        acb_poly_one(ref);
+
+        for (slong i = 0; i < ctx->npts; i++)
+        {
+            acb_poly_zero(ref + 1);
+            acb_poly_set_coeff_acb(ref + 1, 0, ctx->pts + i);
+            acb_poly_set_coeff_si(ref + 1, 1, 1);
+            acb_poly_atan_series(ref + 1, ref + 1, ctx->nder, ctx->prec);
+
+            for (slong m = 0; m < dop_order; m++)
+            {
+                if (!acb_poly_overlaps(ref + m, acb_ode_sol_sum_ptr(ctx->sol + m, i, 0)))
+                {
+                    flint_printf("FAIL\n\n");
+                    flint_printf("i = %wd, m = %wd, a = %{acb}\n", i, m, ctx->pts + i);
+                    for (slong j = 0; j < 2; j++)
+                        flint_printf("sum[%wd] = %{acb_poly}\n", j,
+                                     acb_ode_sol_sum_ptr(ctx->sol + m, i, 0));
+                    flint_printf("ref = %{acb_poly}\n", ref + m);
+                    flint_abort();
+                }
+            }
+        }
+
+        if (ctx->flags & acb_ode_WANT_SERIES)
+        {
+            for (int m = 0; m < dop_order; m++)  /* XXX tmp */
+                acb_poly_truncate(ctx->sol[m].series, nterms);
+
+            acb_poly_zero(ref + 1);
+            acb_poly_set_coeff_si(ref + 1, 1, 1);
+            acb_poly_atan_series(ref + 1, ref + 1, nterms, ctx->prec);
+
+            for (slong m = 0; m < dop_order; m++)
+            {
+                if (!acb_poly_overlaps(ref + m, ctx->sol[m].series))
+                {
+                    flint_printf("FAIL\n\n");
+                    flint_printf("m = %wd, nterms = %wd\n", m, nterms);
+                    flint_printf("series = %{acb_poly}\n", ctx->sol[m].series);
+                    flint_printf("ref = %{acb_poly}\n", ref);
+                    flint_abort();
+                }
+            }
+
+        }
+
+        _acb_poly_vec_clear(ref, dop_order);
+
+        acb_ode_sum_context_clear(ctx);
+    }
+
+    TEST_FUNCTION_END(state);
+}
diff --git a/src/acb_poly.h b/src/acb_poly.h
index eb2c8565ee..dcfa1de04a 100644
--- a/src/acb_poly.h
+++ b/src/acb_poly.h
@@ -684,6 +684,20 @@ acb_poly_allocated_bytes(const acb_poly_t x)
     return _acb_vec_allocated_bytes(x->coeffs, x->alloc);
 }
 
+/* Vector functions */
+
+acb_poly_struct * _acb_poly_vec_init(slong n);
+void _acb_poly_vec_clear(acb_poly_struct *vec, slong n);
+void _acb_poly_vec_zero(acb_poly_struct *dest, slong n);
+void _acb_poly_vec_set(acb_poly_struct *dest, const acb_poly_struct *src, slong n);
+void _acb_poly_vec_set_block(acb_poly_struct *dest, const acb_poly_struct *src,
+                             slong n, slong base, slong len);
+void _acb_poly_vec_fit_length(acb_poly_struct *vec, slong n, slong len);
+void _acb_poly_vec_set_length(acb_poly_struct *vec, slong n, slong len);
+void _acb_poly_vec_normalise(acb_poly_struct *vec, slong n);
+int _acb_poly_vec_overlaps(acb_poly_struct *vec1, acb_poly_struct *vec2, slong n);
+
+
 #ifdef __cplusplus
 }
 #endif
diff --git a/src/acb_poly/vec.c b/src/acb_poly/vec.c
new file mode 100644
index 0000000000..25658ce525
--- /dev/null
+++ b/src/acb_poly/vec.c
@@ -0,0 +1,80 @@
+#include "acb_poly.h"
+
+acb_poly_struct *
+_acb_poly_vec_init(slong n)
+{
+    acb_poly_struct *vec = (acb_poly_struct *) flint_malloc(sizeof(acb_poly_struct) * n);
+
+    for (slong i = 0; i < n; i++)
+        acb_poly_init(vec + i);
+
+    return vec;
+}
+
+void
+_acb_poly_vec_clear(acb_poly_struct *vec, slong n)
+{
+    for (slong i = 0; i < n; i++)
+        acb_poly_clear(vec + i);
+    flint_free(vec);
+}
+
+void
+_acb_poly_vec_zero(acb_poly_struct *dest, slong n)
+{
+    for (slong i = 0; i < n; i++)
+        acb_poly_zero(dest + i);
+}
+
+void
+_acb_poly_vec_set(acb_poly_struct *dest, const acb_poly_struct *src, slong n)
+{
+    for (slong i = 0; i < n; i++)
+        acb_poly_set(dest + i, src + i);
+}
+
+void
+_acb_poly_vec_set_block(acb_poly_struct *dest, const acb_poly_struct *src,
+                        slong n, slong start, slong len)
+{
+    for (slong k = 0; k < n; k++)
+    {
+        if (start >= (src + k)->length)
+            continue;
+        slong len1 = FLINT_MIN(len, (src + k)->length - start);
+        acb_poly_fit_length(dest + k, len1);
+        _acb_vec_set((dest + k)->coeffs, (src + k)->coeffs + start, len1);
+        _acb_poly_set_length(dest + k, len1);
+        _acb_poly_normalise(dest + k);
+    }
+}
+
+void
+_acb_poly_vec_fit_length(acb_poly_struct *vec, slong n, slong len)
+{
+    for (slong i = 0; i < n; i++)
+        acb_poly_fit_length(vec + i, len);
+}
+
+void
+_acb_poly_vec_set_length(acb_poly_struct *vec, slong n, slong len)
+{
+    for (slong i = 0; i < n; i++)
+        _acb_poly_set_length(vec + i, len);
+}
+
+void
+_acb_poly_vec_normalise(acb_poly_struct *vec, slong n)
+{
+    for (slong i = 0; i < n; i++)
+        _acb_poly_normalise(vec + i);
+}
+
+int
+_acb_poly_vec_overlaps(acb_poly_struct *vec1, acb_poly_struct *vec2, slong n)
+{
+    for (slong i = 0; i < n; i++)
+        if (!acb_poly_overlaps(vec1 + i, vec2 + i))
+            return 0;
+    return 1;
+}