# Nicer way to plug an R function into the RHS of an ODEProblem?
**URL:** https://discourse.julialang.org/t/nicer-way-to-plug-an-r-function-into-the-rhs-of-an-odeproblem/114640
**Category:** Modelling & Simulations
**Created:** [May 23, 2024, 7:22pm UTC](https://discourse.julialang.org/t/nicer-way-to-plug-an-r-function-into-the-rhs-of-an-odeproblem/114640 "2024-05-23T19:22:43Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@sdwfrost](https://discourse.julialang.org/u/sdwfrost)
#### Post date: [May 23, 2024, 7:22pm UTC](https://discourse.julialang.org/t/nicer-way-to-plug-an-r-function-into-the-rhs-of-an-odeproblem/114640/1 "2024-05-23T19:22:43Z")
Hi All,
I’m trying to plug an R function into an ODEProblem in OrdinaryDiffEq. The below code works, but is there a simple way to make the code neater e.g. define ODEProblem in terms of an RObject?
```Julia
using OrdinaryDiffEq
using RCall
using Plots
R"""
sir_ode_r <- function(u,p,t){
S <- u[1]
I <- u[2]
R <- u[3]
N <- S+I+R
beta <- p[1]
cee <- p[2]
gamma <- p[3]
dS <- -beta*cee*I/N*S
dI <- beta*cee*I/N*S - gamma*I
dR <- gamma*I
return(c(dS,dI,dR))
}
"""
function sir_ode_jl(u,p,t)
robj = rcall(:sir_ode_r, u, p, t)
return convert(Array,robj)
end
δt = 0.1
tmax = 40.0
tspan = (0.0,tmax)
u0 = [990.0,10.0,0.0] # S,I,R
p = [0.05,10.0,0.25] # β,c,γ
prob = ODEProblem(sir_ode_jl, u0, tspan, p)
sol = solve(prob, Tsit5(), dt = δt)
plot(sol)
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 23, 2024, 7:26pm UTC](https://discourse.julialang.org/t/nicer-way-to-plug-an-r-function-into-the-rhs-of-an-odeproblem/114640/2 "2024-05-23T19:26:28Z")
For this kind of function, you can use ModelingToolkit to translate it to a symbolic form and generate the Julia function. That would also fix the performance. That’s how the diffeqr/diffeqpy JIT compiler works.
Otherwise, I think this is how I would expect it to look.
---
### Author:  [@sdwfrost](https://discourse.julialang.org/u/sdwfrost)
#### Post date: [May 23, 2024, 7:43pm UTC](https://discourse.julialang.org/t/nicer-way-to-plug-an-r-function-into-the-rhs-of-an-odeproblem/114640/3 "2024-05-23T19:43:12Z")
Cool - can you show a simple example (modeltoolkitize?)
Would this also work with my PythonCall example?
```julia
using OrdinaryDiffEq
using PythonCall
using Plots
@pyexec """
def sir_ode_py(u,p,t):
S = u[0]
I = u[1]
R = u[2]
N = S+I+R
beta = p[0]
c = p[1]
gamma = p[2]
dS = -beta*c*I/N*S
dI = beta*c*I/N*S - gamma*I
dR = gamma*I
return [dS,dI,dR]
""" => sir_ode_py
sir_ode_jl(u,p,t) = pyconvert(Array{Float64}, sir_ode_py(u, p, t))
δt = 0.1
tmax = 40.0
tspan = (0.0,tmax)
u0 = [990.0,10.0,0.0] # S,I,R
p = [0.05,10.0,0.25] # β,c,γ
prob = ODEProblem{false}(sir_ode_jl, u0, tspan, p)
sol = solve(prob, Tsit5(), dt = δt)
plot(sol)
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 23, 2024, 7:44pm UTC](https://discourse.julialang.org/t/nicer-way-to-plug-an-r-function-into-the-rhs-of-an-odeproblem/114640/4 "2024-05-23T19:44:57Z")
> [@sdwfrost](#):
>
> Would this also work with my PythonCall example?
>
> ```julia
>
> ```
It should. Try just calling `modelingtoolkitize(prob)`. If you make that convert into an `Array` instead of forcing `Array{Float64}` you should be good.
---
### Author:  [@sdwfrost](https://discourse.julialang.org/u/sdwfrost)
#### Post date: [May 23, 2024, 7:58pm UTC](https://discourse.julialang.org/t/nicer-way-to-plug-an-r-function-into-the-rhs-of-an-odeproblem/114640/5 "2024-05-23T19:58:09Z")
For Python, this code works (converting to `Array` rather than `Array{Float64}`
```julia
using ModelingToolkit
@named sys = modelingtoolkitize(prob)
prob_mtk = ODEProblem(sys, u0, tspan, p)
sol_mtk = solve(prob_mtk, Tsit5(), dt = δt)
plot(sol_mtk)
```
However, modelingtoolkitize breaks for the R version: `ERROR: MethodError: no method matching sexpclass(::Num)`
(yes, I know I’m not using MTK9 🙂 )
---
### Author:  [@sdwfrost](https://discourse.julialang.org/u/sdwfrost)
#### Post date: [May 23, 2024, 8:15pm UTC](https://discourse.julialang.org/t/nicer-way-to-plug-an-r-function-into-the-rhs-of-an-odeproblem/114640/6 "2024-05-23T20:15:44Z")
PS. Is there any way to pass the names of states and parameters when using `modelingtoolkitize` so that the resulting equations are easier to read?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [May 23, 2024, 9:26pm UTC](https://discourse.julialang.org/t/nicer-way-to-plug-an-r-function-into-the-rhs-of-an-odeproblem/114640/7 "2024-05-23T21:26:59Z")
No, but it wouldn’t take more than 15 minutes to add. Open an issue so hopefully I don’t forget. At a conference right now but it’ll end up in the email list that way.