# Solving for the steady-state of a large nonlinear ODE in SciML

**URL:** <https://discourse.julialang.org/t/solving-for-the-steady-state-of-a-large-nonlinear-ode-in-sciml/50613>\
**Category:** Modelling & Simulations\
**Tags:** sciml, nlsolve, newton, nonlinear, steady-state\
**Created:** [November 23, 2020, 5:59am UTC](https://discourse.julialang.org/t/solving-for-the-steady-state-of-a-large-nonlinear-ode-in-sciml/50613 "2020-11-23T05:59:21Z")\
**Posts on this page:** 20\
**Page:** 1

<div class="post-metadata">

**Author:** ![briochemc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/briochemc/32/4209_2.png) [@briochemc](https://discourse.julialang.org/u/briochemc)\
**Post date:** [November 23, 2020, 5:59am UTC](https://discourse.julialang.org/t/solving-for-the-steady-state-of-a-large-nonlinear-ode-in-sciml/50613/1 "2020-11-23T05:59:21Z")

</div>

I have a few questions:

1. What is the difference in scope between NonlinearSolve.jl and NLSolvers.jl/NLsolve.jl?
2. What are the pros/cons of NonlinearSolve.jl vs NLSolvers.jl/NLsolve.jl?
3. Are there any plans to include the algorithms of SIAMFANLEquations.jl into either/both of NonlinearSolve.jl and NLSolvers.jl/NLsolve.jl?
4. In particular, any plans to include the pseudo-transient continuation solver? Or maybe there are some other versions out there?
5. Are there other packages/solvers I should know about for solving large nonlinear systems?

For context, I deal with large-ish (up to ~1M) systems of nonlinear ODEs of marine tracers, for which I develop AIBECS.jl to, well, simulate these tracers. AIBECS essentially provides the tooling for users to create the functions that define the system of ODEs. In SciML’s lingo, it creates two functions of the state variable `u` and the parameters `p`:

```julia
F(u,p) = f(u, p, t) -> the rate of change of the system
∇ₓF(u,p) = jac(u, p, t) -> the jacobian of f with respect to u

```

As you may notice, there is no `t` dependency for my case of `f` and `jac`, and so most of the time, I am looking for the steady-state of these ODEs. To solve for that steady-state, I currently use my own rewrite of @ctkelley’s MATLAB `nsold.m` (which is probably not optimal, in the sense that my rewrite is not optimal, not his algorithms, to be precise! 😅). In other words, for a given set of parameters `p`, I feed `x -> f(x,p,0)` and `x -> jac(x,p,0)` to the newton solver to find the root `x` such that `f(x,p,0) ≈ 0`. Now this all works already quite well in my opinion, but I’m always looking to improve the package, and at the time of writing, it seems I have many choices to replace my non-expert code, among which

- [NLsolve.jl](https://github.com/JuliaNLSolvers/NLsolve.jl) by @pkofod, the most used package for solving large nonlinear systems in Julia AFAIK
- [NLSolvers.jl](https://github.com/JuliaNLSolvers/NLSolvers.jl), @pkofod’s own rewrite of NLsolve.jl, destined to replace it, but already available for those that want to live on the edge
- JuliaComputing/SciML’s [NonlinearSolve.jl](https://github.com/JuliaComputing/NonlinearSolve.jl), which seems to have a large overlap with NLsolve.jl/NLSolvers.jl although it seems to be targeting specifically SciML problems
- [SIAMFANLEquations.jl](https://github.com/ctkelley/SIAMFANLEquations.jl), which is @ctkelley’s own Julia version of `nsold.m` and co, being developed in parallel to writing the corresponding book. Particularly of interest to me is the pseudo-transient continuation solver, which I don’t think has any equivalent in other Julia packages, and is very relevant to my research, where the initial guess may be quite far and impede the Newton algorithm convergence.

Please forgive me if I butchered your work with my ignorant definitions and assumptions, I have the utmost respect for the developers/scientists involved with these and I think these are all amazing pieces of work/software. But back to my questions: I’m not really sure what I should be using in the future and keep asking myself these or similar questions, in my head or on slack, and I thought I should lay these down here for the experts (@pkofod, @ctkelley, @ChrisRackauckas and all the relevant SciML people, and whoever I should know about) to chime in if they want! 🙂

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [November 23, 2020, 7:58am UTC](https://discourse.julialang.org/t/solving-for-the-steady-state-of-a-large-nonlinear-ode-in-sciml/50613/2 "2020-11-23T07:58:05Z")

</div>

You might be interested in [BifurcationKit.jl](https://github.com/rveltz/BifurcationKit.jl) which is geared towards large scale problems. You have a Newton solver which works with sparse matrix or krylov methods for any precision (Float32…) and has been tested on GPU. Additionally, you have continuation routines which allow you to use homotopy methods.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 23, 2020, 8:38am UTC](https://discourse.julialang.org/t/solving-for-the-steady-state-of-a-large-nonlinear-ode-in-sciml/50613/3 "2020-11-23T08:38:32Z")

</div>

SciML made a common interface for LinearProblem, NonlinearProblem, QuadratureProblem, and OptimizationProblem. The reason is because this is required in order to be a build target for ModelingToolkit, and also to allow for swapping the libraries internally (as you’ve seen with SteadyStateDiffEq.jl!). In the place of Quadrature.jl is a universal wrapper over the other quadrature libraries, which also adds AD support. GalacticOptim.jl does the same for OptimizationProblem. NonlinearSolve.jl is meant to do the same.

The integration with ModelingToolkit makes a lot of sense because that means you can use symbolic tooling to enhance the numerical speeds. [https://github.com/JuliaComputing/StructuralTransformations.jl](https://github.com/JuliaComputing/StructuralTransformations.jl) already has nonlinear tearing for MTK-specific `NonlinearSystem`s, so if you go through MTK you can end up with a much simpler nonlinear solve. You also get automated parallelism, sparsity detection, analytical Jacobians, etc. There’s more that can be done there too, with BLT transformations and automation of program transformation to enforce some nonlinear constraints.

NonlinearSolve.jl started by implementing a bunch of scalar nonlinear problems. The reason is because… Roots.jl wasn’t worth wrapping due to its allocations even when scalar. But next up we’re going to have this library `@requires` wrap NLsolve, NLsolvers, SUNDIALS KINSOL, SIAMFANLEquations.jl, etc. and everything else we can find. Then we add AD support. Then we integrate it into tooling like SteadyStateDiffEq.jl.

The reason is because then you can just pass `NLsolveAlg()` into a bigger algorithm and it’ll know how to dispatch the internal nonlinear solver. This will make algorithms built on nonlinear solvers (the problem we run into) be a lot more flexible without the extra work on the user. So, the same kind of deal as we’ve done with differential equations (and optimization and quadrature).

But we’ll leave the implementation of nonlinear solvers to the ecosystem. That’s @pkofod and @ctkelley’s job, so we’ll make use of those. I don’t think we’ll get into the business of making special nonlinear solvers… I think 😉.

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [November 23, 2020, 11:41am UTC](https://discourse.julialang.org/t/solving-for-the-steady-state-of-a-large-nonlinear-ode-in-sciml/50613/4 "2020-11-23T11:41:04Z")

</div>

I will be registering [SIAMFANLEquations.jl](https://github.com/ctkelley/SIAMFANLEquations.jl) soon. The only thing left to do is get the first part of Chapter 2 into the [notebook](https://github.com/ctkelley/NotebookSIAMFANL) . That part will have an example of how to use `ptcsol.jl` (the pseudo-transient continuation code) for the buckling beam problem.

`ptcsol.jl` and `nsol.jl` (nsold.m with some new things like mixed-precision solves) are ready to go. @briochemc, I’d be happy to assist you in your application of `ptcsol.jl`. If you are already using it, please tell me how it’s going.

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [November 23, 2020, 12:38pm UTC](https://discourse.julialang.org/t/solving-for-the-steady-state-of-a-large-nonlinear-ode-in-sciml/50613/5 "2020-11-23T12:38:11Z")

</div>

@ctkelley Sorry for my ad, I just note the reference of my package from yours.

@ChrisRackauckas One additional point is whether we want to support types that are not `<: AsbtractArray`

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [November 23, 2020, 12:55pm UTC](https://discourse.julialang.org/t/solving-for-the-steady-state-of-a-large-nonlinear-ode-in-sciml/50613/6 "2020-11-23T12:55:28Z")

</div>

> [@rveltz](#):
>
> @ChrisRackauckas One additional point is whether we want to support types that are not `<: AsbtractArray`

Yes, if the solvers do. And they already do: `Number` is not an `AbstractArray`.

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [November 23, 2020, 1:09pm UTC](https://discourse.julialang.org/t/solving-for-the-steady-state-of-a-large-nonlinear-ode-in-sciml/50613/7 "2020-11-23T13:09:26Z")

</div>

@rveltz, no need to apologize. The more data we as a collective can give @briochemcm, the better off he will be.

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [November 23, 2020, 1:10pm UTC](https://discourse.julialang.org/t/solving-for-the-steady-state-of-a-large-nonlinear-ode-in-sciml/50613/8 "2020-11-23T13:10:20Z")

</div>

Yup, scalar equation solvers would get very upset if fed an array. @rveltz, do you have something in mind?

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [November 23, 2020, 2:05pm UTC](https://discourse.julialang.org/t/solving-for-the-steady-state-of-a-large-nonlinear-ode-in-sciml/50613/9 "2020-11-23T14:05:41Z")

</div>

Things like ArrayFire, ApproxFun…

---

<div class="post-metadata">

**Author:** ![briochemc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/briochemc/32/4209_2.png) [@briochemc](https://discourse.julialang.org/u/briochemc)\
**Post date:** [November 30, 2020, 6:21am UTC](https://discourse.julialang.org/t/solving-for-the-steady-state-of-a-large-nonlinear-ode-in-sciml/50613/10 "2020-11-30T06:21:23Z")

</div>

Thank you! After looking a bit more carefully at the code, I realize that there is no option for an out-of-place Jacobian function, so I need to spend a bit more time with the inner workings of AIBECS to return an in-place Jacobian function before I can try `ptcsol` (forming the Jacobian out-of-place was never the bottleneck for me so I just skipped trying to implement that until now).

---

<div class="post-metadata">

**Author:** ![briochemc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/briochemc/32/4209_2.png) [@briochemc](https://discourse.julialang.org/u/briochemc)\
**Post date:** [November 30, 2020, 8:00am UTC](https://discourse.julialang.org/t/solving-for-the-steady-state-of-a-large-nonlinear-ode-in-sciml/50613/11 "2020-11-30T08:00:53Z")

</div>

Actually, I realize now that updating a sparse Jacobian inplace is not so straightforward… See, e.g. these many discussions around this topic:

- [Update matrix block inplace from a sparse array](https://discourse.julialang.org/t/update-matrix-block-inplace-from-a-sparse-array/36360),
- [Efficiently updating values in a sparse diagonal matrix](https://discourse.julialang.org/t/efficiently-updating-values-in-a-sparse-diagonal-matrix/13111),
- [Update bands of a CSC Matrix](https://discourse.julialang.org/t/update-bands-of-a-csc-matrix/36800),
- [Optimal way to do in-place modification of sparse matrix?](https://discourse.julialang.org/t/optimal-way-to-do-in-place-modification-of-sparse-matrix/5170),
- [A version of sparse! that updates a previously existing SparseMatrixCSC?](https://discourse.julialang.org/t/a-version-of-sparse-that-updates-a-previously-existing-sparsematrixcsc/31943)
- [Example use of `sparse!` - #2 by bennedich](https://discourse.julialang.org/t/example-use-of-sparse/21808/2)

* * *

@ctkelley, Are there any plans to allow for out-of-place (sparse) jacobian functions? I ask because I already have available a fast-but-out-of-place jacobian function 🤷‍♂️

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [November 30, 2020, 8:14am UTC](https://discourse.julialang.org/t/solving-for-the-steady-state-of-a-large-nonlinear-ode-in-sciml/50613/12 "2020-11-30T08:14:04Z")

</div>

Hi,

I think you forgot these alternatives:

- use of BandedArrays (can really be much faster than sparse and ~ easy to update)
- use of [SparseDiffTools.jl](https://github.com/JuliaDiff/SparseDiffTools.jl) which can be inplace

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [November 30, 2020, 8:25am UTC](https://discourse.julialang.org/t/solving-for-the-steady-state-of-a-large-nonlinear-ode-in-sciml/50613/13 "2020-11-30T08:25:19Z")

</div>

> [@briochemc](#):
>
> After looking a bit more carefully at the code, I realize that there is no option for an out-of-place Jacobian function

What do you mean? You can’t cook an inplace one from your out-of place one?

---

<div class="post-metadata">

**Author:** ![mauro3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mauro3/32/292_2.png) [@mauro3](https://discourse.julialang.org/u/mauro3)\
**Post date:** [November 30, 2020, 8:33am UTC](https://discourse.julialang.org/t/solving-for-the-steady-state-of-a-large-nonlinear-ode-in-sciml/50613/14 "2020-11-30T08:33:03Z")

</div>

First create a sparse matrix with “filled-locations” wherever the Jacobian can be non-zero. Then during your in-place update, fill those locations. If one of them is ==0, that’s ok, just set it to 0 (and make sure that the method you use for setting it, does not delete the “filled-location” spot).

But yes, consider SparseDiffTools.jl

---

<div class="post-metadata">

**Author:** ![ctkelley](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ctkelley/32/10684_2.png) [@ctkelley](https://discourse.julialang.org/u/ctkelley)\
**Post date:** [November 30, 2020, 11:22am UTC](https://discourse.julialang.org/t/solving-for-the-steady-state-of-a-large-nonlinear-ode-in-sciml/50613/15 "2020-11-30T11:22:26Z")

</div>

Once you’ve allocated the storage, ptcsol defaults to the usual sparse factorization from SparseSuite. I’ll have an example of this in a later chapter of the notebook.

For now I’m using special structrues. There’s an example for PTC and the buckling beam in

> **[GitHub - ctkelley/NotebookSIAMFANL: Notebook for my new solver book](https://github.com/ctkelley/NotebookSIAMFANL)**
>
> Notebook for my new solver book. Contribute to ctkelley/NotebookSIAMFANL development by creating an account on GitHub.

where the Jacobian is tridiagonal. There is also an Newton’s method example using BandedMatrices.jl. You need to know the sparsity pattern in advance so the solver does not allocate for the Jacobian. Let me think a bit about how hard it would be for me to use an allocating Jacobian evaluation function.

@rveltz was correct on both of his points. If you Jacobian is banded, using something like BandedMatrices.jl will make your life better and sparse differencing can figure out the sparsity pattern automatically.

I will tag the latest version of the notebook today and will register the package by the end of this week.

---

<div class="post-metadata">

**Author:** ![briochemc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/briochemc/32/4209_2.png) [@briochemc](https://discourse.julialang.org/u/briochemc)\
**Post date:** [November 30, 2020, 12:52pm UTC](https://discourse.julialang.org/t/solving-for-the-steady-state-of-a-large-nonlinear-ode-in-sciml/50613/16 "2020-11-30T12:52:38Z")

</div>

Well I think I could, it’s just a bit more complicated than I anticipated because my out-of-place Jacobian function drops zeros (the non-structural entries that are sometimes zero). This is both good (because there are quite a few) and bad (because some branching inside the nonlinear function means some non-structural zeros appear and disappear and make it hard to use for an inplace version).

For the record and because it helps me think about, I’ll desribe my Jacobian function here and I might eventually copy this into the package docs anyway. Also maybe someone has good ideas to improve it and make it inplace, too, so here goes. This goes in the details so it is a little wordy — so apologies in advance.

TLDR: I have a function F that looks like F(x) = G(x) - T x where G is a perfect candidate for ForwardDiff and SparseDiffTools, while T is already available as the Jacobian of x \mapsto T x. How can I make an inplace Jacobian for J(x) = \nabla G(x) - T?

* * *

Long version:

When trying to solve for F(x) = 0, the variable x is actually a concatenation of smaller (but still large) vectors:

x = \begin{bmatrix} x\_1\\ x\_2\\ \vdots \\ x\_N \\ \end{bmatrix}

where each x\_i can be a ~200,000 sized column vector.

All the x\_i have the same length and share the same indices, in the sense that these indices correspond to some linear indices of a 3D grid. Each different x\_i thus actually represents a different 3D field that has been “vectorized”. The function F can be similarly split

F(x) = \begin{bmatrix} F\_1(x\_1, \ldots , x\_N) \\ F\_2(x\_1, \ldots , x\_N) \\ \vdots \\ F\_N(x\_1, \ldots , x\_N) \\ \end{bmatrix}.

and each F\_i is further split into a “nonloncal” linear part and a “local” nonlinear part, like so

F\_i(x) = G\_i(x\_1, \ldots , x\_N) - T\_i x\_i

where the T\_i are sparse matrices but not diagonal, and their sparsity structure does not correspond to any of the sparse array packages I have seen (they have 1 to ~20 off-diagonals that can be quite far off the main diagonal and whose distance to the diagonal can vary with each row/column). The G\_i, on the other hand, are nonlinear, but local, in that the j-th entry [G\_i(x\_1, \ldots , x\_N)]\_j only depends on the j-th entries of its x\_i arguments. The AIBECS.jl user supplies the T\_i and the G\_i functions, which I then use to build F and its Jacobian as a function of a single larger vector x.

The jacobian of F can be built efficiently using autodifferentiation of the G\_i only (ForwardDiff.jl in my case). If I denote by \partial\_j G\_i the “derivative” of G\_i with respect to x\_j to be the (square) jacobian of the x\_j \mapsto G\_i(x\_1,\ldots,x\_N) function, which is diagonal because the G\_i are “local” as explained just above, then jacobian of F can be written as

J(x) = \begin{bmatrix} (\partial\_1 G\_1)(x)-T\_1 & (\partial\_2 G\_1)(x) & \cdots & (\partial\_N G\_1)(x)\\ (\partial\_1 G\_2)(x) & (\partial\_2 G\_2)(x) -T\_2 & \cdots & (\partial\_N G\_2)(x)\\ \vdots & & \ddots & \vdots \\ (\partial\_1 G\_N)(x) & (\partial\_2 G\_N)(x) & \cdots & (\partial\_N G\_N)(x) - T\_N\\ \end{bmatrix}

In practice, I actually build it with

J(x) = \begin{bmatrix} (\partial\_1 G\_1)(x) & (\partial\_2 G\_1)(x) & \cdots & (\partial\_N G\_1)(x)\\ (\partial\_1 G\_2)(x) & (\partial\_2 G\_2)(x) & \cdots & (\partial\_N G\_2)(x)\\ \vdots & & \ddots & \vdots \\ (\partial\_1 G\_N)(x) & (\partial\_2 G\_N)(x) & \cdots & (\partial\_N G\_N)(x)\\ \end{bmatrix} - \begin{bmatrix} T\_1 & & & \\ & T\_2 & & \\ & & \ddots & \\ & & & T\_N\\ \end{bmatrix}

where N ForwardDiff’d calls to the G\_i's are sufficient to build the first term (which is essentially an array of diagonals) with some combination of `sparse`, `Diagonal`, `hcat`, and `vcat` operations, and where the second term uses the `blockdiag` function on the T\_i. I believe the `-` operation drops all the non-structural zeros (but I need to double check), which is good and bad, as mentioned in the intro to this post.

I could use SparseDiffTools for the first term by supplying the sparsity pattern, but it is not efficient to use it on the whole F because the T\_i fill in a lot of the sparsity pattern and I actually don’t need to use ForwardDiff on the T\_i x\_i terms since these are linear and the T\_i are already provided.

So how should I go about making an inplace version of this Jacobian?

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [November 30, 2020, 1:47pm UTC](https://discourse.julialang.org/t/solving-for-the-steady-state-of-a-large-nonlinear-ode-in-sciml/50613/17 "2020-11-30T13:47:59Z")

</div>

Using FD, you can update the non diagonal blocks inplace. So I would focus on updating the block (1,1) inplace. You save the indices of T\_1 which are changed by \partial\_1G\_1 and write them inplace. This is what I do [here](https://github.com/rveltz/BifurcationKit.jl/blob/master/src/periodicorbit/PeriodicOrbitFD.jl#L732)

---

<div class="post-metadata">

**Author:** ![briochemc](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/briochemc/32/4209_2.png) [@briochemc](https://discourse.julialang.org/u/briochemc)\
**Post date:** [February 8, 2021, 5:43am UTC](https://discourse.julialang.org/t/solving-for-the-steady-state-of-a-large-nonlinear-ode-in-sciml/50613/18 "2021-02-08T05:43:03Z")

</div>

@rveltz I just started thinking about all this again and noticed your link does not work anymore. I’m assuming the code was just moved somewhere else — could you update said link? 🙏

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [February 8, 2021, 6:53am UTC](https://discourse.julialang.org/t/solving-for-the-steady-state-of-a-large-nonlinear-ode-in-sciml/50613/19 "2021-02-08T06:53:15Z")

</div>

> [@rveltz](#):
>
> here

id say [https://github.com/rveltz/BifurcationKit.jl/blob/master/src/Utils.jl#L103](https://github.com/rveltz/BifurcationKit.jl/blob/master/src/Utils.jl#L103)

---

<div class="post-metadata">

**Author:** ![rveltz](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rveltz/32/2707_2.png) [@rveltz](https://discourse.julialang.org/u/rveltz)\
**Post date:** [February 8, 2021, 7:05am UTC](https://discourse.julialang.org/t/solving-for-the-steady-state-of-a-large-nonlinear-ode-in-sciml/50613/20 "2021-02-08T07:05:24Z")

</div>

Then it is used here [BifurcationKit.jl/PeriodicOrbitTrapeze.jl at master · bifurcationkit/BifurcationKit.jl · GitHub](https://github.com/rveltz/BifurcationKit.jl/blob/master/src/periodicorbit/PeriodicOrbitTrapeze.jl#L883)

[Next page](https://discourse.julialang.org/t/solving-for-the-steady-state-of-a-large-nonlinear-ode-in-sciml/50613.md?page=2)