import numpy as np
# pip install numpy ,system 
print( "..............................................................................................")
print( " simple reaction :n*y -> products,dy/dt = -k*(y)**n , rk45 ,version 1.0 ")
print( " peter.vlasschaert@gmail.com,16/01/2023                                                       ")
print( "..............................................................................................")
def rk45(f, y0, t0, t_end, h):
    """
    Solves the ODE y'(t) = f(t, y(t)) with initial condition y(t0) = y0 using the RK45 method.    
    """
    t = t0
    y = y0
    solutions = [y] # all output during simulation
    while t < t_end:
        k1 = h * f(t, y)
        k2 = h * f(t + h/4, y + k1/4)
        k3 = h * f(t + 3*h/8, y + 3*k1/32 + 9*k2/32)
        k4 = h * f(t + 12*h/13, y + 1932*k1/2197 - 7200*k2/2197 + 7296*k3/2197)
        k5 = h * f(t + h, y + 439*k1/216 - 8*k2 + 3680*k3/513 - 845*k4/4104)
        k6 = h * f(t + h/2, y - 8*k1/27 + 2*k2 - 3544*k3/2565 + 1859*k4/4104 - 11*k5/40)
        
        y_next = y + 25*k1/216 + 1408*k3/2565 + 2197*k4/4104 - k5/5
        t += h
        y = y_next
        solutions.append(y) # all output during simulation
    return solutions

y0 = 1 # initial concentration ,see figure
t0 = 0
t_end = 1
h = 0.1
k = 1 # rate reaction constant
print("--------------------------------")
print("reaction model : n*y-> products ")
print("dy/dt =-k*y^n,n=integer,y0 @t=0 ")
print("................................")
def create_ode_fn(n, k):
    def ode_fn(t, y):
        return -k*y**n
    return ode_fn


import matplotlib.pyplot as plt
# pip install matplotlib ,system
for n in [0,1,2,3]:
    ode_fn = create_ode_fn(n,k)
    solutions = rk45(ode_fn, y0, t0, t_end, h)
    print("dim : solutions = ",n,len(solutions))
    # Create a list of time points
    times = np.arange(t0, t_end+2*h, h)
    print("dim : times = ",n,len(times))
    # Plot the solutions
    plt.plot(times, solutions, label='n = {}'.format(n))
print(" end : simulation > plot with n values ")
# Add labels and title to graphics
plt.xlabel('Time')
plt.ylabel('Solution')
plt.title('RK45 solutions for different n')
plt.legend()
plt.grid()
# Show the plot output
plt.show()