Homotopy Continuation Method to Find All Root of a Polynomial Equation I
Author's Data:  Housam BINOUS
Department of Chemical Engineering
National Institute of Applied Sciences and Technology
Tunis, TUNISIA

Reference: M. Kuno and J. D. Seader, Computing All Real Solutions to Systems of Nonlinear Equations with a Global Fixed-Point Homotopy,  Ind. Eng. Chem. Res., 1988, 27, 1320-1329.

Solve give the roots of the polynomial equation

(x - 0.25) (x - 0.75) (x - 0.11) (x - 2) // 
      Expand

Solve[
      0.04125` - 0.615625` x + 2.5175` (x^2) - 
          3.11` (x^3) + (x^4) == 0, x]

Fixed Point Homotopy : 
H(x,t)=(1-t) (x-x0) + t f(x) with x0=2.5

H[x_, t_] := (1 - t) (x - 2.5) + 
        t (0.04125` - 0.615625` x + 2.5175` (x^2) - 
              3.11` (x^3) + (x^4))

H[x, t]

Off[General::"spell"]

Grad[fon_, var_List] := D[fon, #] & /@ var

J = Grad[H[x, t], {x, t}]

Js = J /. x -> x[s] /. t -> t[s]

Eq[1] = D[x[s], {s, 1}] ==   First[Drop[Js, 1]]

Eq[2] = D[t[s], {s, 1}] ==  -First[Drop[Js, {2}]]

sol = 
      NDSolve[{Eq[1], Eq[2], t[0] ==  0, x[0] ==  2.5}, {x[s], 
          t[s]}, {s, -4.38, 0.094}, 
        MaxStepFraction ->  0.0005]

As expected we find four solutions at t=1

plt1 = 
      ParametricPlot[{t[s], x[s]} /. sol, {s, -4.3378, 0.093}, 
        Epilog ->  {RGBColor[1, 0, 0], Line[{{1, -3}, {1, 3}}]}, 
        PlotRange -> {{0.98, 1.02}, {-0.5, 2.2}}]

l = 
      Graphics[{{RGBColor[1, 0, 1], PointSize[0.03], 
            Point[{1, 0.25}]}, {RGBColor[1, 0, 1], PointSize[0.03], 
            Point[{1, 0.75}]}, {RGBColor[1, 0, 1], PointSize[0.03], 
            Point[{1, 2}]}, {RGBColor[1, 0, 1], PointSize[0.03], 
            Point[{1, 0.11}]}}]

Show[plt1, l, Frame -> True]

sol1 = 
      FindRoot[(t[s] /. sol /. s ->  a) == 1, {a, 
          0.01}]

Evaluate[x[s] /. sol /. s -> a /. sol1 ]

sol1 = 
      FindRoot[(t[s] /. sol /. s ->  a) == 
          1, {a, -3}]

Evaluate[x[s] /. sol /. s -> a /. sol1 ]

sol1 = 
      FindRoot[(t[s] /. sol /. s ->  a) == 
          1, {a, -3.8, -4}]

Evaluate[x[s] /. sol /. s -> a /. sol1 ]

sol1 = 
      FindRoot[(t[s] /. sol /. s ->  a) == 
          1, {a, -4.3, -4}]

Evaluate[x[s] /. sol /. s -> a /. sol1 ]

We have found all solutions on the same single branch because there are no real roots for the equation: f(x) - x + (Subscript[x, O]) = 0 and x -> ± ∞. This is a result of our choice of (Subscript[x, O]).

Solve[
      0.04125` - 0.615625` x + 
          2.5175` (x^2) - 
3.11` (x^3) + (x^4) - x + 
          2.5 == 0, x]

ParametricPlot[{t[s], x[s]} /. sol, {s, -4.337875792747977`, 
        0.09376277513774145`}, 
      PlotRange -> {{-0.05, 0.02}, {-200, 200}}, 
      PlotPoints ->  10000, 
      PlotStyle ->  RGBColor[1, 0, 0]]

Newton Homotopy :
H(x,t)= f(x) - (1-t) f(x0) with x0= - 0.5

0.04125` - 0.615625` x + 2.5175` (x^2) - 
        3.11` (x^3) + (x^4) /. x ->   -0.5

H[x_, t_] := -(1 - t) (1.4296875) + 
        0.04125` - 0.615625` x + 2.5175` (x^2) - 
        3.11` (x^3) + (x^4)

J = Grad[H[x, t], {x, t}]

Js = J /. x -> x[s] /. t -> t[s]

Eq[1] = D[x[s], {s, 1}] ==   First[Drop[Js, 1]]

Eq[2] = D[t[s], {s, 1}] == - First[Drop[Js, {2}]]

sol = 
      NDSolve[{Eq[1], Eq[2], t[0] ==  0, 
          x[0] ==  -0.5}, {x[s], t[s]}, {s, -0.8, 2}, 
        MaxStepFraction ->  0.001]

As expected we find four solutions at t=1

plt2 = 
      ParametricPlot[{t[s], x[s]} /. sol, {s, -0.2, 1.8}, 
        Epilog ->  {RGBColor[1, 0, 0], Line[{{1, -3}, {1, 3}}]}, 
        PlotRange -> {{0.98, 1.02}, {-0.5, 2.5}}]

l = 
      Graphics[{{RGBColor[1, 0, 1], PointSize[0.03], 
            Point[{1, 0.25}]}, {RGBColor[1, 0, 1], PointSize[0.03], 
            Point[{1, 0.75}]}, {RGBColor[1, 0, 1], PointSize[0.03], 
            Point[{1, 2}]}, {RGBColor[1, 0, 1], PointSize[0.03], 
            Point[{1, 0.11}]}}]

Show[plt1, l, Frame -> True]

sol1 = 
      FindRoot[(t[s] /. sol /. s ->  a) == 1, {a, 
          0.01}]

Evaluate[x[s] /. sol /. s -> a /. sol1 ]

sol1 = 
      FindRoot[(t[s] /. sol /. s ->  a) == 1, {a, 
          0.5}]

Evaluate[x[s] /. sol /. s -> a /. sol1 ]

sol1 = 
      FindRoot[(t[s] /. sol /. s ->  a) == 1, {a, 
          0.9}]

Evaluate[x[s] /. sol /. s -> a /. sol1 ]

sol1 = 
      FindRoot[(t[s] /. sol /. s ->  a) == 1, {a, 1.5, 
          1.8}]

Evaluate[x[s] /. sol /. s -> a /. sol1 ]