i_did_your_homework_for_you

import matplotlib.pyplot as plt
import numpy as np
from scipy.optimize import fsolve


# Functions
def f1(x):
    return x**3 + (c - a) * x**2 / b - 1


def f2(x):
    return (1 - y / y0)**(eta - 1) * x**7 + b * x**6 - a * x**4 - b

# Initialization
k = 0.46235
r0 = 0.42
n0 = 101400 * 4 * np.pi * r0**3 / (3 * 8.3143 * 303)
b = 16 * np.pi * k * r0**2 / (3 * n0)
y0 = 49000
eta = 0.0289 * 9.98 * y0 / (8.3143 * 303)
n = 45
a = n / n0
M = 1.12
c = M / (0.0289 * n0)

lamba = (fsolve(f1, y0)**0.5)[0]
y_max = y0 - ((b + a * lamba**4 - b * lamba**6) / lamba**7)**(1 / (eta - 1)) * y0
print("Max Altitude: {:4.1f}".format(y_max))

h = np.zeros(100)
v = np.zeros(100)
i = 0

for y in np.linspace(0, y_max, 100):
    h[i] = y

    lambb = fsolve(f2, 49000)[0]
    v[i] = ((20 * 9.8 * r0 / 3) * (lambb - c * (1 - y / y0)**(1 - eta) / lambb**2))**0.5
    i += 1

plt.plot(h / 1000, v)
plt.xlabel("Distance x (km)")
plt.ylabel("Velocity v (m/s)")
plt.title("Plot of velocity over distance")
plt.show()