Let’s solve the linear ODE `u'=1.01u`

. First setup the package:

Define the derivative function `f(u,p,t)`

.

Then we give it an initial condition and a time span to solve over:

With those pieces we call `diffeqr::ode.solve`

to solve the ODE:

This gives back a solution object for which `sol$t`

are the time points and `sol$u`

are the values. We can check it by plotting the solution:

Now let’s solve the Lorenz equations. In this case, our initial condition is a vector and our derivative functions takes in the vector to return a vector (note: arbitrary dimensional arrays are allowed). We would define this as:

```
f <- function(u,p,t) {
du1 = p[1]*(u[2]-u[1])
du2 = u[1]*(p[2]-u[3]) - u[2]
du3 = u[1]*u[2] - p[3]*u[3]
return(c(du1,du2,du3))
}
```

Here we utilized the parameter array `p`

. Thus we use `diffeqr::ode.solve`

like before, but also pass in parameters this time:

```
u0 = c(1.0,0.0,0.0)
tspan <- list(0.0,100.0)
p = c(10.0,28.0,8/3)
sol = diffeqr::ode.solve(f,u0,tspan,p=p)
```

The returned solution is like before. It is convenient to turn it into a data.frame:

Now we can use `matplot`

to plot the timeseries together:

Now we can use the Plotly package to draw a phase plot:

Plotly is much prettier!

If we want to have a more accurate solution, we can send `abstol`

and `reltol`

. Defaults are `1e-6`

and `1e-3`

respectively. Generally you can think of the digits of accuracy as related to 1 plus the exponent of the relative tolerance, so the default is two digits of accuracy. Absolute tolernace is the accuracy near 0.

In addition, we may want to choose to save at more time points. We do this by giving an array of values to save at as `saveat`

. Together, this looks like:

```
abstol = 1e-8
reltol = 1e-8
saveat = 0:10000/100
sol = diffeqr::ode.solve(f,u0,tspan,p=p,abstol=abstol,reltol=reltol,saveat=saveat)
udf = as.data.frame(sol$u)
plotly::plot_ly(udf, x = ~V1, y = ~V2, z = ~V3, type = 'scatter3d', mode = 'lines')
```