I know there is a “convert_system” function which can convert a “NonlinearSystem” to an “ODESystem”. Can we do it backwards? So that we can have both steady-state and dynamic simulation with a single code base.

Good point. Open an issue. Right now, if you have a NonlinearSystem you can use the SteadyStateDiffEq.jl methods to solve the nonlinear system via ODEs, but it’s not at the symbolic level.

Okay cool. Yeah we’ll get to it: it’s not too difficult. Also the NonlinearSolve.jl docs are what’s on my table right now, and there’s some simple fixes there that would make this straight forward. So check back there in a little bit, and remind me if I drop it.

I am also interested in converting an ODESystem to a NonlinearSystem. I wonder if this feature has been implemented or if there is a good workaround @ChrisRackauckas. I am aware of doing something like this for an ODESystem but then we don’t have the nice solution interface anymore

```
dfg = NonlinearProblem(odeprob)
nsol = solve(dfg, SciMLNLSolve.NLSolveJL())
```

The solution interface is the same? MWE?

Here is an MWE showcasing what would be convenient but doesn’t work currently. It boils down to that for an `ODEProblem`

, you can access the solution with `sol[u(t,x,y)]`

while for a `NonlinearProblem`

created from that `ODEProblem`

you cannot access it via `nsol[u(t,x,y)]`

or `nsol[u(x,y)]`

. I could rewrite the equation to a nonlinear problem which would allow to access `nsol[u(x,y)]`

, but it would be nice if there was a convenience function that would do that conversion. I believe that is also what the original post was looking for.

```
using ModelingToolkit, OrdinaryDiffEq, MethodOfLines, Plots, SciMLNLSolve
@parameters x y t
@variables u(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
Dt = Differential(t)
# Boundary conditions
bcs = [
u(0,x,y) ~ 0.0,
u(t,0,y) ~ 0.0, u(t,1,y) ~ 0.0,
u(t,x,0) ~ 0.0, u(t,x,1) ~ 0.0
]
# Space and time domains
domains = [
t ∈ IntervalDomain(0.0,1.0),
x ∈ IntervalDomain(0.0,1.0),
y ∈ IntervalDomain(0.0,1.0)
]
# Discretization
dx = 0.04
nd = round(Int, 1 / dx) + 1
discretization = MOLFiniteDifference([x=>dx,y=>dx],t)
eq = Dt(u(t,x,y)) ~ Dxx(u(t,x,y)) + Dyy(u(t,x,y)) + sin(pi*x)*sin(pi*y)
@named pde_system = PDESystem([eq],bcs,domains,[t,x,y],[u(t,x,y)])
prob = discretize(pde_system,discretization)
sol = solve(prob,TRBDF2(),saveat=0.1)
# solution interface accessible via u(t,x,y)
sol[u(t,x,y)]
dfg = NonlinearProblem(prob)
nsol = solve(dfg, SciMLNLSolve.NLSolveJL())
# solution interface not accessible via u(t,x,y): ArgumentError: u(t, x, y) is neither an observed nor a state variable.
nsol[u(t,x,y)] # does not work
# instead, we have to access the solution via the solution vector - this is not convenient for larger models
M_ss = zeros(nd,nd)
M_ss[2:end-1, 2:end-1] .= reshape(nsol.u, nd-2, nd-2)
```

MWE adapted from Control of PDE systems · JuliaSimControl (juliahub.com)

@ChrisRackauckas - you can find the MWE in the post above. And btw, big thanks to all your great work with SciML and MethodOfLines.jl in particular

Yeah that’s Unable to Index into a steadystate solution using symbols · Issue #338 · SciML/MethodOfLines.jl · GitHub, it’s on our mind and needs to get fixed.