# Edit:

So the way I proposed doing it here *does* work, I just made a mistake that @Oscar_Smith pointed out, so using the functions this way is fine.

# Original Q:

Hi there!

I’m still a bit new to Julia and trying to figure out the best way of using `DifferentialEquations.jl`

to set up a large system of ODEs for simulation. I have a large system of equations (over 50) that represent biomass flux between containers and I’m looking for a way to define that system via a set of functions, because each container is technically defined by the same equation, just with different varying parameters.

As an example with only two containers, imagine a different system:

I know that I can set that up in Julia with the `DifferentialEquations.jl`

package like this:

```
r1 = 1.5
r2 = 1.9
a12 = 0.010975
a21 = 0.01
m1 = 1.4
m2 = 1.3
k1 = 5000
k2 = 2000
function eqns!(dx, x, p, t)
dx[1] = r1*x[1]*(1-(x[1]/k1)) + x[1]*(a12*x[2])*(1-(x[1]/k1)) -
x[1]*(a21*x[2]) - m1*x[1]
dx[2] = r2*x[2]*(1-(x[2]/k2)) + x[2]*(a21*x[2])*(1-(x[2]/k2)) -
x[2]*a12*x[1] - m2*x[2]
end
x0 = [100.0; 1000.0]
tspan = (0.0, 10.0)
prob = ODEProblem(eqns!, x0, tspan)
sol = solve(prob)
plot(sol)
```

Which gives the resulting plot:

But in this example I have to define every parameter manually and then write out each equation. I want a way to set each parameter programatically so that in a 50-container system, rather than manually setting a_{ij} manually for each combo of *i* and *j* where i \text{ and } j = \{x \space \epsilon \space \mathbb{Z} \space | \space 0 < x \leq 50 \} I would simply call some function `get_a(i, j)`

and get that value for any *i* and *j* combo.

I’m thinking of this in the context of turning each term into a function so that for any X_i, the eq’n would be:

\frac{dX_i}{dt} =

`indep_gain(i)`

+ `dep_gain(i)`

- `dep_removal(i)`

- `indep_removal(i)`

For my example, I thought of doing that this way where I define a series of functions (this is just a small play example that I came up with) and then the functions can be called for any X_i:

```
function indep_growth(i)
#define r
if i == 1
r = 1.5
k = 5000
else
r = 1.8
k = 2000
end
# return product of r and abundance of focal x
return (r*x[i](1-(x[i]/k)))
end
function dep_growth(i)
# determine if a12 or a21 should be used
if i == 1 # use a12
a = 010975
k = 5000
return x[i]*(a*x[2])*(1-(x[i]/k))
else
a = 0.01
k = 2000
return x[i]*(a*x[1])*(1-(x[i]/k))
end
end
function dep_death(i)
if i == 1
death = dep_growth(2)
else
death = dep_growth(1)
end
return death
end
function indep_death(i)
if i == 1
m = 1.4
else
m = 1.3
end
return m*x[i]
end
function eqns!(dx, x, p, t)
for i in 1:2
dx[i] = indep_growth(i) + dep_growth(i) - dep_death(i) - dep_death(i)
end
end
```

And then creating the `eqns!()`

with something like a loop:

```
function eqns!(dx, x, p, t)
for i in 1:2
dx[i] = indep_growth(i) + dep_growth(i) - dep_death(i) - dep_death(i)
end
end
```

But however, when I then try to use `solve(prob)`

, I get

```
ERROR: MethodError: objects of type Float64 are not callable
```

With a huge stacktrace.

**So my question is**: Is there a way to assign the actual equations in the function via a series of functions as opposed to writing them all out manually? I have re-generate my 50-eqn system like 1000 times so I can’t write it out manually each time.

Any thoughts/help is much appreciated!!