Solve quintic equation using differential equation

Question

Background

DiffResolvent.nb, here gives a method that solve quintic equations using differential equations.

There's a progress to transform polynomial equation into differential equation.

Then he gave an example:

$$x^5-5 x^3+5 x-t = 0$$

The roots of the equation satisfy the differential equation:

$$25 \left(4-t^2\right) x''(t)-25 t x'(t)+x(t)=0$$

Solve that and then can get the roots of quintic equation.

$$x(t)=c_1 \cos \left(\frac{1}{5} \arcsin\frac{t}{2}\right)-c_2 \sin \left(\frac{1}{5} \arcsin\frac{t}{2}\right)$$

We can easily check the results:

eq=x^5-5x^3+5 x-t
resolvent=25 (4-t^2) x''[t]-25 t x'[t]+x[t]

dso=DSolveValue[resolvent==0,x[t],t]/.Sqrt[4-t^2]->I Sqrt[t^2-4]
ser=Series[Evaluate[eq/.x->dso],{t,0,1},Assumptions->t\[Element]Reals];
sol=dso/.Solve[CoefficientList[ser,t]==0,{C[1],C[2]}];

SortBy[x/.NSolve[eq/.t->3,x],Im]
SortBy[sol/.t->3//N//Chop,Im]

I think this can't be solved by Solve.

Root[-3+5 #1-5 #1^3+#1^5&,1]//ToRadicals

---

Problem

But these codes were written in Mathematica 3.0 and can't run any more.

The behaviour of these functions differs from DifferentialRoot.

There's no intermediate process so I can't find whats wrong with these codes.

Goal

A function which convert the polynomial equation into differential equation.

Answer 1

First of all, the code doesn't even work in v3:

After checking the notebook further, I found the code line

algeqn = Collect[Numerator[algeqn], t[rho] ] 

quite suspicious, because

1. It doesn't seem to match the corresponding description

   We replace powers t[rho]^k with exponents k greater than 4 using the original quintic equation

2. It actually doesn't do anything because algeqn is an equation but Numerator doesn't have effect on equation (at least from v3).

I believe the algorithm should be coded as e.g.:

eqn = t[rho]^5 - t[rho] - rho == 0

diffeqn = a1 t''''[rho] + a2 t'''[rho] + a3 t''[rho] + 
    a4 t'[rho] + a5 t[rho] + a6  == 0

deriv = Flatten[Table[Solve[D[eqn, {rho, k}], D[t[rho], {rho, k}]], {k, 1, 4}]]

algeqn = Simplify[diffeqn //. deriv]
(* Here's the key point: *)
expr = FixedPoint[
  Collect[#, t@rho] /. t[rho]^i_ /; i > 4 :> (t[rho] + rho) t[rho]^(i - 5) &, 
  Numerator@Together[Subtract @@ algeqn]]

coe = Solve[
  CoefficientList[expr, t[rho]] == 0 // Thread, {a1, a2, a3, a4, a5, a6}]

sol = First@DSolve[diffeqn /. coe // Simplify, t@rho, rho]; // AbsoluteTiming
(* {70.096278, Null} *)
approximation = sol /. HoldPattern@ HypergeometricPFQ[w__] -> 1

eqnapprox = eqn /. approximation

system = (#1 == 0 &) /@ Take[CoefficientList[eqnapprox[[1]], rho], 4]

coeC = Solve[system, C /@ Range@4]

solfinal = sol /. coeC
(* Check: *)
Block[{rho = RandomReal[1, WorkingPrecision -> 16]}, eqn /. solfinal]