#225 Finite difference method to solve ODE using python

Let y(t) be the altitude of the missile at time t. The gravity g = 9.8 m/s^2. We want to have the missile at 50 m off the ground after t = 5 s after launch.
Find the velocity at launch v0=y'(0)? Ignore the drag of the air resistance.
Boundary conditions are y(0) = 0 and y(5) = 50.
d^2y/dt^2=-g
The result should be y'(0)=34.5 m/s.

d^2y/dt^2=(y(i-1) – 2*y(i) + y(i+1))/(h^2)

Since the time interval is [0,5] and we have n = 100, h = (5-0)/n, using the finite difference approximated derivatives, we obtain
y(0) = 0
y(i−1) − 2y(i) + y(i+1) = −g*h^2, (i = 1,2,…,n − 1),
y(n) = 50

import numpy as np
import matplotlib.pyplot as plt
plt.style.use("seaborn-poster")
n = 100
h = (5-0) / n
# Get A
A = np.zeros((n+1, n+1))
A[0, 0] = 1
A[n, n] = 1

for i in range(1, n):
    A[i, i-1] = 1
    A[i, i] = -2
    A[i, i+1] = 1
print(A)

# Get b
b = np.zeros(n+1)
b[1:-1] = -9.8*h**2
b[-1] = 50
print(b)

# solve the linear equations
y = np.linalg.solve(A, b)
t = np.linspace(0, 5, n+1)
plt.figure(figsize=(10,8))
plt.plot(t, y)
plt.plot(5, 50, "ro")
plt.xlabel("time (s)")
plt.ylabel("altitude (m)")
plt.savefig('ex225.png', dpi=72)
plt.show()
y_n_1 = -9.8*h**2 + 2*y[0] - y[1]
v0=(y[1] - y_n_1) / (2*h)
print(v0)