I’m surprised that I actually found something scientific-computing related that is easier for me to do in JavaScript than it is in Julia. I’ve written functions that use the 4/5th order Runge-Kutta-Fehlberg method to solve a variety of different ODE systems. I’ve been trying to write a more general version of such a function (general in the sense that it should be able to solve ODE systems of any size), namely using the code:

```
using LinearAlgebra;
struct solObj
t
vars
end
function f(t, vars)
x = vars[1];
y = vars[2];
array = [x*cos(y), cos(x)];
return array;
end
function rkf45(f::Function, t0::Number, tf::Number, conds, epsilon::Float64, dtInitial::Float64)
dt = dtInitial;
t = Float64[t0];
vars = conds;
i = 1;
while ( t[i] < tf )
dt = min(dt, tf-t[i]);
K1 = dt*f(t[i], vars[i,:]);
K2 = dt*f(t[i] + dt/4, vars[i,:] + K1/4);
K3 = dt*f(t[i] + 3*dt/8, vars[i,:] + 3*K1/32 + 9*K2/32);
K4 = dt*f(t[i] + 12*dt/13, vars[i,:] + 1932*K1/2197 - 7200*K2/2197 + 7296*K3/2197);
K5 = dt*f(t[i] + dt, vars[i,:] + 439*K1/216 - 8*K2 + 3680*K3/513 - 845*K4/4104);
K6 = dt*f(t[i] + dt/2, vars[i,:] - 8*K1/27 + 2*K2 - 3544*K3/2565 + 1859*K4/4104 - 11*K5/40);
vars1 = vars[i,:] + 25*K1/216 + 1408*K3/2565 + 2197*K4/4104 - K5/5;
vars2 = vars[i,:] + 16*K1/135 + 6656*K3/12825 + 28561*K4/56430 - 9*K5/50 + 2*K6/55;
R = min(abs.(vars2-vars1)/dt...);
s = (epsilon/(2*R))^(1/4);
if ( R <= epsilon )
push!(t, t[i]+dt);
push!(vars, vars1);
i += 1;
end
dt *= s;
end
solution = solObj(t, vars);
return solution;
end
```

(in this case I’ve written a function for `f`

both to give you an example of what my input function could look like and also to help me test out my function). Now, K1, K2, K3, K4, K5 and K6 should ideally all be row vectors, as should vars[i,:], vars1 and vars2. vars is an array with each row corresponding to a different time value and the columns corresponding to different dependent variables (in this case x and y). Vars is used instead of x and y directly because this function is designed to be more generalized, such that it should be able to work on an ODE system of any size.

This code looked to me like it should work fine, but I get the error:

```
ERROR: DimensionMismatch("dimensions must match: a has dims (Base.OneTo(1),), b has dims (Base.OneTo(2),), mismatch at 1")
Stacktrace:
[1] promote_shape at ./indices.jl:178 [inlined]
[2] promote_shape at ./indices.jl:169 [inlined]
[3] +(::Array{Float64,1}, ::Array{Float64,1}) at ./arraymath.jl:45
[4] rkf45(::typeof(f), ::Float64, ::Float64, ::Array{Float64,1}, ::Float64, ::Float64) at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/RKF45.jl:21
[5] top-level scope at REPL[114]:1
[6] run_repl(::REPL.AbstractREPL, ::Any) at /build/julia/src/julia-1.5.2/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:288
```

when I execute the command `rkf45(f, 0.0, 10.0, [0, 0], 1e-6, 0.1)`

. Line 21 is the one on which K1 is defined. Despite this, given the error I’m getting I’m guessing it’s actually an issue with line 22 (where K2 is defined), as there’s actually an addition there. I’ve tried resolving the issue using the `reshape()`

function, but feels like playing whac-a-mole as whenever one issue is temporarily resolved with it another pops up. For example the code:

```
using LinearAlgebra;
struct solObj
t
vars
end
function f(t::Number, vars)
x = vars[1];
y = vars[2];
array = [x*cos(y), cos(x)];
return array;
end
function rkf45(f::Function, t0::Number, tf::Number, conds::Array, epsilon::Float64, dtInitial::Float64)
dt = dtInitial;
t = [t0];
N = length(conds);
vars = reshape(conds, 1, N);
i = 1;
while ( t[i] < tf )
dt = min(dt, tf-t[i]);
K1 = reshape(dt*f(t[i], reshape(vars[i,:], 1, N)), 1, N);
K2 = reshape(dt*f(t[i] + dt/4, reshape(vars[i,:], 1, N) + K1/4), 1, N);
K3 = reshape(dt*f(t[i] + 3*dt/8, reshape(vars[i,:], 1, N) + 3*K1/32 + 9*K2/32), 1, N);
K4 = reshape(dt*f(t[i] + 12*dt/13, reshape(vars[i,:], 1, N) + 1932*K1/2197 - 7200*K2/2197 + 7296*K3/2197), 1, N);
K5 = reshape(dt*f(t[i] + dt, reshape(vars[i,:], 1, N) + 439*K1/216 - 8*K2 + 3680*K3/513 - 845*K4/4104), 1, N);
K6 = reshape(dt*f(t[i] + dt/2, reshape(vars[i,:], 1, N) - 8*K1/27 + 2*K2 - 3544*K3/2565 + 1859*K4/4104 - 11*K5/40), 1, N);
vars1 = reshape(vars[i,:], 1, N) + 25*K1/216 + 1408*K3/2565 + 2197*K4/4104 - K5/5;
vars2 = reshape(vars[i,:], 1, N) + 16*K1/135 + 6656*K3/12825 + 28561*K4/56430 - 9*K5/50 + 2*K6/55;
R = min(abs.(vars2-vars1)/dt...);
s = (epsilon/(2*R))^(1/4);
if ( R <= epsilon )
push!(t, t[i]+dt);
push!(vars, vars1);
i += 1;
end
dt *= s;
end
solution = solObj(t, vars);
return solution;
end
```

when I run `rkf45(f, 0.0, 10.0, [0, 0], 1e-6, 0.1)`

yields:

```
ERROR: MethodError: no method matching push!(::Array{Int64,2}, ::Array{Float64,2})
You might have used a 2d row vector where a 1d column vector was required.
Note the difference between 1d column vector [1,2,3] and 2d row vector [1 2 3].
You can convert to a column vector with the vec() function.
Closest candidates are:
push!(::Any, ::Any, ::Any) at abstractarray.jl:2252
push!(::Any, ::Any, ::Any, ::Any...) at abstractarray.jl:2253
push!(::Base.Nowhere, ::Any) at array.jl:1347
...
Stacktrace:
[1] rkf45(::typeof(f), ::Float64, ::Float64, ::Array{Int64,1}, ::Float64, ::Float64) at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/RKF45-reshape.jl:34
[2] top-level scope at REPL[6]:1
[3] run_repl(::REPL.AbstractREPL, ::Any) at /build/julia/src/julia-1.5.2/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:288
```

I have tried using `push!(vars, vec(vars1)')`

instead like this error suggests may be best (as instead of a column vector, I want it to be a row vector) but that just gives the error:

```
ERROR: MethodError: no method matching push!(::Array{Int64,2}, ::Adjoint{Float64,Array{Float64,1}})
You might have used a 2d row vector where a 1d column vector was required.
Note the difference between 1d column vector [1,2,3] and 2d row vector [1 2 3].
You can convert to a column vector with the vec() function.
Closest candidates are:
push!(::Any, ::Any, ::Any) at abstractarray.jl:2252
push!(::Any, ::Any, ::Any, ::Any...) at abstractarray.jl:2253
push!(::Base.Nowhere, ::Any) at array.jl:1347
...
Stacktrace:
[1] rkf45(::typeof(f), ::Float64, ::Float64, ::Array{Int64,1}, ::Float64, ::Float64) at /data/fusion809/GitHub/mine/scripts/julia-scripts/ODEs/RKF45-reshape.jl:34
[2] top-level scope at REPL[3]:1
[3] run_repl(::REPL.AbstractREPL, ::Any) at /build/julia/src/julia-1.5.2/usr/share/julia/stdlib/v1.5/REPL/src/REPL.jl:288
```

Any ideas of a way around this conundrum?