First Order Differential Equation Solving Algorithm

First order differential equations problems can be written as the following initial value problem:

• x' = f(t, x), 0 < t <= T
• x0 = x(0) [initial value]

Let x[i] = x(i * dt) for i = 0 to i = N - 1, where dt = T / (N - 1). For dt 'small' enough, we can approximate the derivative x'(t) using the forward difference:

x'(t) = (x(t + dt) - x(t)) / dt = f(t, x)
x(t + dt) = x(t) + f(t, x) * dt

Then, we can compute all values of x[i] with the following:

x[i] = { x[i - 1] + f(i * dt, x) * dt    0 < i <= N - 1
       { x0                              i = 0

So, the algorithm simplifies to the following:

x = np.zeros(N)
x[0] = x0
for i in range(1, N):
  x[i] = x[i - 1] + f((i - 1) * dt, x[i - 1]) * dt

[full code implementation]