diff --git a/homotopy.jl b/homotopy.jl
new file mode 100644
index 0000000000000000000000000000000000000000..a5710e918d8ceb60e6b5b06812a6f6e06c5cb492
--- /dev/null
+++ b/homotopy.jl
@@ -0,0 +1,142 @@
+# Path following,
+using LinearAlgebra
+using DifferentialEquations
+using ForwardDiff
+using Plots
+using Printf
+
+function Path(F)
+    """
+    path = Path(F)
+    Build the path following 
+    
+    Input:
+    F : Homotopy function
+        val = F(x,λ)
+        Input:
+            x : Float(n)
+            λ : homotopic parameter Float
+        Output:
+            val : value of the Homotopy, Float(n)
+    Output:
+        H : pathfollowing function
+        sol = H(x0, λ0, λf; sf=1e4, abstol=1e-8, reltol=1e-8, display=true)
+            Input:
+                x0 : initial point, Float(n)
+                λ0 : initial value of the homotopic parameter, Float
+                λf : final value of the homotopic parameter, Float
+
+            Optional input
+                sf : max value of sf, Float
+                abstol : Absolute tolerance for the numerical integration, Float
+                reltol : Relative tolerance for the numerical integration, Float
+                display : print result at each integration, Boolean
+    
+            Output:
+                sol : path, solution of the (IVP)
+    """       
+    
+    jac(f, x) = ForwardDiff.jacobian(f, x)
+    
+    function H(x0, λ0, λf; sf=1e4, abstol=1e-8, reltol=1e-8, display=true, saveat=[], ftol_proj=1e-8)
+        n = size(x0, 1)
+        first = true
+        #old_dc = Array{typeof(x0[1])}(undef, n+1)
+        old_dc = zeros(n+1)
+
+        function rhs!(dc, c, par, s)
+        """
+           Compute de right hand side of the IVP
+           rhs!(dc, c, par, s)
+           Input
+               c : state, Float(n)
+               par : parameter (λ0, λf)
+               s : independent parameter
+           Output
+               dc : tangent vector, Float(n)   
+        """
+            λ0, λf = par
+            n     = size(c, 1)-1
+            g(c)  = F(c[1:n], c[n+1])
+            dF    = jac(g, c)      # dF is the Jacobian Matrix of F
+
+            Q, R  = qr(dF')        # QR decomposation
+
+            dc[:] = Q[:,n+1]       # dc[:] is a vector of norm 1 of the kernel
+            # \dot{\lambda} must be of the same signe than (\lambda_f-\lambda0)) the fisrt time
+            if first                                
+                dc[:] = sign((λf-λ0)*dc[end])*dc
+                first = false
+            # test the direction of the tangent vector for taking the good one
+            else
+                dc[:] = sign(dot(old_dc,dc))*dc   
+            end
+            old_dc[:] = dc
+        
+        end
+    
+        function rhsu!(du, u, p, t)
+            rhs!(du, u, p, t)
+        end
+        
+        # ode problem
+        c0  = [x0; λ0]
+        ode = ODEProblem{true}(rhsu!, c0, (0., sf), [λ0; λf])
+        
+        # callback: termination
+        condition(c,t,integrator) = c[end]-λf
+        affect!(integrator) = terminate!(integrator)
+        cbt  = ContinuousCallback(condition,affect!)
+
+        # callback: projector
+        function g!(out,c,p,t)
+            out[1:n] = F(c[1:n], c[n+1])
+        end
+        cbp = ManifoldProjection(g!, nlopts=Dict([(:ftol, ftol_proj)]))
+        
+        # callback print
+        iter = 1
+        function cb_display(c, t, integrator)
+            @printf("%10d", iter)
+            @printf("%16.8e", norm(F(c[1:n], c[n+1])))
+            @printf("%16.8e", norm(c[1:n]))
+            @printf("%16.8e", c[n+1])
+            println()
+            iter = iter + 1
+        end
+        cbd = FunctionCallingCallback(cb_display)
+        
+        #
+        if display
+            cb = CallbackSet(cbt, cbp, cbd)
+            
+            # init print
+            println("\n     Calls  |F(x,λ)|        |x|             λ \n")
+        
+        else
+            cb = CallbackSet(cbt, cbp)
+        end
+        
+        # resolution
+        if saveat==[]
+            sol = solve(ode, Tsit5(), abstol=abstol, reltol=reltol, callback=cb)
+        else
+            sol = solve(ode, Tsit5(), abstol=abstol, reltol=reltol, callback=cb, saveat=saveat)
+        end
+        
+        if display
+            cf = sol[:,end]
+            print("\n Results of the path solver method:\n")
+            println(" xf         = ", cf[1:n]);
+            println(" λf         = ", cf[1+n]);
+            println(" |F(xf,λf)| = ", norm(F(cf[1:n], cf[1+n])));
+            println(" steps      = ", size(sol.t, 1));
+            println(" status     = ", sol.retcode);
+        end
+        
+        return sol
+    end
+        
+    return H
+  
+end
\ No newline at end of file
diff --git a/transfert-conso-min.mp4 b/transfert-conso-min.mp4
new file mode 100644
index 0000000000000000000000000000000000000000..ba49b749aa324d5dedd6ceca98aff343a6c0d680
Binary files /dev/null and b/transfert-conso-min.mp4 differ