Case B

In this situation, we proceed to eliminate all shifts in so as to obtain a purely differential relation.

> eq[1]:=op(select(has,sys[2],diff(f(n,u),u$3)));

> eq[2]:=op(select(has,sys[2],D[2](f)(n+1,u)));

> eq[3]:=op(select(has,sys[2],f(n+2,u)));

[Procedure to replace var in expr with its value solved from relation.]

> plug:=proc(relation,expr,var)
eval(subs(var=solve(relation,var),expr))
end:

[Procedure that collects in each function call.]

> normalize:=proc(expr)
collect(expr,indets(expr,function),distributed,factor)
end:

> eq[4]:=normalize(plug(diff(eq[1],u),eq[2],D[2](f)(n+1,u)));

> eq[5]:=normalize(plug(eq[1],eq[4],f(n+1,u)));

The following is a purely differential equation satisfied by the integral.

> eq[5]:=normalize(primpart(subs(f(n,u)=f(u),eq[5]),diff(f(u),u$4)));

There only remains to solve it. First, a general solution is obtained by a call to Maple's general dsolve.

> dsol:=op(2,dsolve(eq[5],f(u)));

Next, we use initial conditions.

> solve({seq(subs(n=i,dsol)=simplify(int(subs(n=i,expr),z=0..infinity),symbolic),i=1..4)},{_C1,_C2,_C3,_C4});

> normal(subs(%,dsol),expanded);