-- Classical RK4
-- Steve Kifowit
-- June 17, 2009

-- User input
function f(x,y)
   return x*y^2
end

x0 = 0
y0 = 1
Xend = 1
steps = 20
-- No user input required beyond this point, --
-- unless you want to edit the format strings --

h = ( Xend - x0 ) / steps
xx, yy = x0, y0
print("Starting values:")
form1 = "When x = %.6f,   y = %.6f"
io.write( string.format( form1, xx, yy ), "\n"  )

file = io.open("rk4.dat", "w+")
file:write( xx, ", ", yy, "\n" )

for i = 1, steps do
   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
   file:write( xx, ", ", yy, "\n" )
end

file:close()

print("\nEnding values:")
form2 = "When x = %.6f,   y is approximately %.6f"
io.write( string.format( form2, xx, yy ), "\n" )
print("\nData from individual steps written to rk4.dat")