# -----------------------------------------------
# Lecture 21 - Example - ODEs and Explicit Euler
#   ChEn 263 -- Numerical Tools
#   Douglas R. Tree
# -----------------------------------------------

import numpy as np
import matplotlib.pyplot as plt

# Solve the ODE IVP (i.e. rate equation):
#   dc/dt = -c/tau
#   c(0) = c_0
#
#   tau = 0.6
#   c_0 = 1
#
# Explicit Euler method:
#   c_n+1 = c_n + dt*(-c_n/tau)

# set up time grid
dt = 0.1
t_beg = 0
t_end = 5
nt = int((t_end-t_beg)/dt+1)
t = np.zeros(nt)
c = np.zeros(nt)

# constants in our equation
c0 = 1
tau = 0.6

# Loop to code Explicit Euler
t[0] = t_beg
c[0] = c0
for i in range(nt-1):
    c[i+1] = c[i] + dt*(-c[i]/tau)
    t[i+1] = t[i] + dt
      
# exact solution (for comparison)
t_ex = np.linspace(t_beg, t_end, 101)
c_ex = c0 * np.exp(-t_ex/tau)


# plot the result
plt.rc('font', size=20)
plt.figure(figsize=(9,6))
plt.plot(t, c, 'o--', markersize=9, label='Explicit Euler')
plt.plot(t_ex, c_ex, '-', linewidth=3, label='Exact Solution')
plt.xlabel('t')
plt.ylabel('c')
plt.legend()
plt.show()