# Modified Euler's 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 ):
    w = f(xx,yy)
    yt = yy + h * w
    xt = xx + h
    yy = yy + h * ( w + f(xt,yt) ) / 2
    xx = xt
    xy.append([xx,yy])

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

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