Derive a holonomic system for the inner double integral  of the expection of an Euler characteristic number

<< RISC`HolonomicFunctions`; (* load Christoph's package *)

HolonomicFunctions Package version 1.7.3 (21-Mar-2017)
written by Christoph Koutschan
Copyright Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

ann1 = <<ec1_ann1.m; (* input a D-finite system of the inner single integral *)

ann1 = ann1 /. {s1->2,s2->1,m11->1,m21->-1,m22->1}; (* specify the parameters in ann1 as that in Example 1 *)

supp1 = {1, Der[b], Der[d], Der[b]^2, Der[b] Der[d], Der[d]^2, Der[d]^3} (* a support set for an annihilator of the inner double integral *)

Timing[op1=FindCreativeTelescoping[ann1,Der[s], Support -> supp1][[1,1]];] (* compute an annihilator op1 of the inner double integral with supports in supp1 *)

supp2 = {1, Der[b], Der[d], Der[b]^2, Der[b] Der[d], Der[d]^2, Der[b]^3, Der[b]^2 Der[d], Der[b] Der[d]^2, Der[d]^3};(* another support set for another annihilator of the inner double integral *)

Timing[op2=FindCreativeTelescoping[ann1,Der[s], Support -> supp2][[1,1]];](* compute another annihilator op2 of the inner double integral with supports in supp2 *)

Timing[ann2 = OreGroebnerBasis[{op1, op2}];] (* compute a D-finite system ann2 of the inner double integral *)

Length[UnderTheStaircase[ann2]] (* ann2 has holonomic rank 6*)

ann2 >> ec1_ann2.m; (* output ann2 *)