# Solve 2nd order differential equation using Euler method
# Python 3 version
import math
import numpy as np
#  equation of motion
def dvdt(t,x):
    return -4.0*x
#  analytic solution
def sol(t,x0):
    return x0 * math.cos (2.0 * t)
# set parameters
t0  =0.0        # starting time
tmax=10.0       # ending time
dt  =0.01       # delta t
nend=int(tmax/dt)+1
outputstep = 10 # output data every xx steps
# data for plotting graph
i=0
iend=int(tmax/dt/outputstep)+1
t_plot = np.zeros(iend)
x_plot = np.zeros(iend)
a_plot = np.zeros(iend)
#
x0=1.0   # initial value x(0)
v0=0.0   # initial value v(0)
# initial set up
t=t0
x=x0
dxdt=v0
t_plot[i] = t
x_plot[i] = x
a_plot[i] = sol(t,x0)
# first line of the data
print('{:^12}{:^12}{:^12}{:^12}'.format('i   t','numerical','analytic','diff'))
print(f'{i:4d}{t:6.2f}{x:12.5f}{sol(t,x0):12.5f}{x-sol(t,x0):12.5f}')
# t-loop
for n in range(1, nend):
#   forward difference
    dxdt = dxdt + dt * dvdt(t,x)
    x    = x    + dt * dxdt
#   next t
    t = t + dt
    if np.mod(n,outputstep) == 0:
        i=i+1
        t_plot[i] = t
        x_plot[i] = x
        a_plot[i] = sol(t,x0)
        print(f'{i:4d}{t:6.2f}{x:12.5f}{sol(t,x0):12.5f}{x-sol(t,x0):12.5f}')




import matplotlib.pyplot as plt
fig=plt.subplots(figsize=(8,6))
plt.plot(t_plot, x_plot,'.', c='k',label='numerical',markersize=6)
plt.plot(t_plot, a_plot, c='r',label='analytic')
plt.xlabel('t')
plt.ylabel('x')
plt.grid()
plt.legend(loc='lower right')