function f(x)
  return [x[2], (1. - x[1]^2) * x[2] - x[1]]
end


function ode(; x0=[1., 0.], dt=1e-3, t_max=10., dim=2)
  ts = 0. : dt : t_max
  xs = zeros(dim, size(ts)[1])

  x = x0
  xs[:, 1] = x

  for i in 2 : size(ts)[1]
    dx = f(x) .* dt
    x += dx
    xs[:, i] = x
  end
  return (ts, xs)
end


function main()
  t_max = 10.
  @time ts, xs = ode(t_max=t_max)
end


main()