MAXDEORDER=5

Clear[HolonomicDE];
HolonomicDE[f_,F_[x_],print_:False]:=Module[{DF,eqcl,rat,A,approach,sol,eq,i,j,l,m},
DF=Table[D[f,{x,i}],{i,0,MAXDEORDER-1}];
DF=Map[If[Head[#]===Plus,Apply[List,#],{#}]&,DF];(* DF consists of lists of derivatives (sums are lists of summands) *)
eqcl={A[1]DF[[1]]};
For[i=2,i≤MAXDEORDER,i++,
For[l=1,l≤Length[DF[[i]]],l++, (* for every summand of the (i-1)-th derivative do the following *)
m=Length[eqcl];
For[j=1,j≤m,j++,
rat=Together[DF[[i,l]]/eqcl[[j,1]]]; (* building quotient of l-th summand and first element of j-th eq. class*)
If [PolynomialQ[Denominator[rat],x] && PolynomialQ[Numerator[rat],x], (* test, if rat is rational or not  *)
AppendTo[eqcl[[j]],A[i]DF[[i,l]]]; (* add to equivalence class *)
Break[],
If[j==Length[eqcl],AppendTo[eqcl,{A[i]DF[[i,l]]}];  (* build new equivalence class *)
];
];
];
];
var=Table[A[k],{k,i}];
approach=Map[Apply[Plus,#]&,eqcl]; (* sum over all elements in equivalence classes *)
approach=Map[#==0&,approach];
If[print,Print[{eqcl,approach}]];
Off[Solve::"svars"];
sol=Solve[approach,var];
If[sol=!={{}} && Union[var/.sol[[1]]]=!={0},
eq=Sum[A[k+1]D[F[x],{x,k}],{k,0,i-1}]/.sol[[1]];
eq=Numerator[Together[eq]]/.A[_]->1;
Return[Collect[eq,D[F[x],{x,_}],Factor]==0]];
];
]

Aufgabe 1

Clear[DEtoRE];
DEtoRE[DE_,F_[x_],A_[n_]]:=Module[{expr,rules},expr=Expand[DE[[1]]-DE[[2]]];
rules={x^j_ F[x]->A[n-j],x F[x]->A[n-1],F[x]->A[n],x^j_ D[F[x],{x,k_}]->Pochhammer[n+1-j,k] A[n+k-j],x D[F[x],{x,k_}]->Pochhammer[n,k] A[n+k-1],D[F[x],{x,k_}]->Pochhammer[n+1,k] A[n+k]};
Collect[expr/.rules,A[_],Factor]==0]

(a)

DEtoRE[D[f[x],{x,2}]-x f[x]==0,f[x],A[k]]

(b)

DEtoRE[(1-x^2)D[f[x],{x,2}]-2x D[f[x],x]+n^2f[x]==0,f[x],A[k]]

HolonomicDE

MAXDEORDER=5

Clear[HolonomicDE];
HolonomicDE[f_,F_[x_],print_:False]:=Module[{DF,eqcl,rat,A,approach,sol,eq,i,j,l,m},
DF=Table[D[f,{x,i}],{i,0,MAXDEORDER-1}];
DF=Map[If[Head[#]===Plus,Apply[List,#],{#}]&,DF];(* DF consists of lists of derivatives (sums are lists of summands) *)
eqcl={A[1]DF[[1]]};
For[i=2,i≤MAXDEORDER,i++,
For[l=1,l≤Length[DF[[i]]],l++, (* for every summand of the (i-1)-th derivative do the following *)
m=Length[eqcl];
For[j=1,j≤m,j++,
rat=Together[DF[[i,l]]/eqcl[[j,1]]]; (* building quotient of l-th summand and first element of j-th eq. class*)
If [PolynomialQ[Denominator[rat],x] && PolynomialQ[Numerator[rat],x], (* test, if rat is rational or not  *)
AppendTo[eqcl[[j]],A[i]DF[[i,l]]]; (* add to equivalence class *)
Break[],
If[j==Length[eqcl],AppendTo[eqcl,{A[i]DF[[i,l]]}];  (* build new equivalence class *)
];
];
];
];
var=Table[A[k],{k,i}];
approach=Map[Apply[Plus,#]&,eqcl]; (* sum over all elements in equivalence classes *)
approach=Map[#==0&,approach];
If[print,Print[approach]];
Off[Solve::"svars"];
sol=Solve[approach,var];
If[sol=!={{}} && Union[var/.sol[[1]]]=!={0},
eq=Sum[A[k+1]D[F[x],{x,k}],{k,0,i-1}]/.sol[[1]];
eq=Numerator[Together[eq]]/.A[_]->1;
Return[Collect[eq,D[F[x],{x,_}],Factor]==0]];
];
]

(c)

DE=HolonomicDE[Sin[x]+ArcTan[x],f[x]]

DEtoRE[DE,f[x],A[k]]

Aufgabe 2

Prozedur "SumToHyper"

Clear[SumToHyper];
SumToHyper[a_,k_]:=Module[{rat,i,l,m,p,upp,low,arg,coeff,start},
rat=FunctionExpand[(a/.k->k+1)/a];
l=FactorList[rat];
low={};upp={};arg=1;
For[i=1,i≤Length[l],i++,
If[FreeQ[l[[i]],k],
arg*=l[[i,1]]^l[[i,2]],
coeff=Coefficient[l[[i,1]],k,1];
If [l[[i,2]]>0,
Do[AppendTo[upp,Coefficient[l[[i,1]],k,0]/coeff];
arg*=coeff,{l[[i,2]]}],
Do[AppendTo[low,Coefficient[l[[i,1]],k,0]/coeff];
arg/=coeff,{-l[[i,2]]}];
];
];
];
start=0;
m=Min[Select[low,IntegerQ] ];
If[IntegerQ[m],
upp=Map[#-m+1&,upp];
low=Map[#-m+1&,low];
start+=-m+1;
];
p=Min[Select[upp,IntegerQ] ];
If[IntegerQ[p] && p<1,
upp=Map[#-p+1&,upp];
low=Map[#-p+1&,low];
start+=-p+1;
];
If[MemberQ[low,1],
low=Delete[low,Position[low,1][[1]]],
AppendTo[upp,1]
];
Return[Simplify[a/.k->start]Hypergeom[upp,low,arg]];
]

(a)

a=(-1)^k x^(2 k)/(2 k)!

pFq=SumToHyper[a,k]

Apply[HypergeometricPFQ,pFq]

(b)

a=Binomial[n,k] x^k

pFq=SumToHyper[a,k]

Apply[HypergeometricPFQ,pFq]

(c)

a=(-1)^(k+1)/k x^k

pFq=SumToHyper[a,k]

pFq/.Hypergeom[x_,y_,z_]->HypergeometricPFQ[x,y,z]

(d)

a=(-1)^k x^(2 k+1)/(2 k+1)

pFq=SumToHyper[a,k]

pFq/.Hypergeom[x_,y_,z_]->HypergeometricPFQ[x,y,z]