# Classic 4th-order Runge-Kutta method for dy/dx=f(x,y), y(x0)=y0.
# Initial and ending values are shown, as is the scatterplot of approximations.
# To see the entire list of values, uncomment the line near the bottom.
# You can run this code in a SageMath cell.

from math import *
import matplotlib.pyplot as plt

#---------------------------------
def f(x,y):
    return(x+y)
#---------------------------------
x0 = 1.0
y0 = 2.0
Xend = 2.0
steps = 10
#----------------------------------
# No user input required beyond this point

h = ( float(Xend) - float(x0) ) / steps
xx, yy = float(x0), float(y0)
xy = [[xx,yy]]


print("Starting values:")
print("x = ", xx, ", y = ",  yy)

for i in range( steps ):
    k1 = h * f(xx,yy)
    ww = yy + k1 / 2
    k2 = h * f(xx+h/2,ww)
    ww = yy + k2 / 2
    k3 = h * f(xx+h/2,ww)
    ww = yy + k3
    k4 = h * f(xx+h,ww)
    yy = yy + ( k1 + 2*k2 + 2*k3 + k4 ) / 6
    xx = xx + h
    xy.append([xx,yy])

print("\nEnding values:")
print("x = ", xx, ", y = ",  yy)

#print("\n", *xy, sep="\n" )
plt.scatter(*zip(*xy))
plt.show()