We are going to show a roundabout but rather effective method of solving some rather complicated (definite) integrals in closed-form via a rather unusual method: a special factorization of linear ODEs.  The example we will use is a 2 week old question that has yet to get an answer.

First, the problem: compute the integral

f  := exp(a*m)* (sqrt(T)*sqrt(Pi)*erf((1/2)*sqrt(2)*(m+a*T)/sqrt(T))*exp(a*m)*a -sqrt(T)*sqrt(Pi)*exp(a*m)*a +sqrt(2)*exp(-(1/2)*(m^2+a^2*T^2)/T)):  gg := Int(f,m=0..K);

The first step is to translate this into a problem on ODEs:

dd := PDEtools[dpolyform](y(K)=gg,no_Fn);
dd2 := numer((lhs-rhs)(op([1,1],dd)));

which gives a 4th order linear ODE with polynomial coefficients. Then we need to draw on the power of the 'diffop' (for differential operators) sub-package of DEtools, in particular the LCLM factorization.

oper := DEtools[de2diffop](dd2, y(K), [Dx,x]);
fact := DEtools[DFactorLCLM](oper,[Dx,x]);
de3 := map(DEtools[diffop2de], fact, [Dx,x], y(x));

So now we have something quite fantastic: the solution to the integral will be made up of the solutions to those 2 second-order ODEs. What will those look like?

basis4  := map(op,map(dsolve,de3,output=basis));
eq := add(C[i]*basis4[i],i=1..4);

Those are exactly as asked by the original poster! We need to find the exact solution, so we just form the linear system from the first 4 Taylor coefficients, solve it, and substitute

eqs := [seq(eval(diff(eq,[x$i]),x=0) - eval(value(diff(gg,[K$i])),K=0),i=0..3)]:
solve(eqs,{C[1],C[2],C[3],C[4]});
eval(eq,%);