-- Butcher's 5th-order Runge-Kutta method for solving
-- dy/dx = f(x,y), y(x0)=y0

-- Steve Kifowit
-- October 21, 2009

-- f is defined here.  Math functions require "math."
-- For a list of available math functions, see the Lua Reference 
-- Manual available at http://www.lua.org.
---------------------------------
function f(x,y)
   return 1/2 - y/(10+4*x)
end
---------------------------------

x0 = 0        -- Initial x value
y0 = 0        -- Initial y value
Xend = 10     -- Ending x value
steps = 40    -- Number of steps
----------------------------------

-- 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("rk5.dat", "w+")
file:write( xx, ", ", yy, "\n" )

for i = 1, steps do
   k1 = h * f(xx,yy)
   ww = yy + k1 / 4
   k2 = h * f(xx+h/4,ww)
   ww = yy + (k1 + k2) / 8
   k3 = h * f(xx+h/4,ww)
   ww = yy + k3 - k2 / 2
   k4 = h * f(xx+h/2,ww)
   ww = yy + (3*k1 + 9*k4) / 16
   k5 = h * f(xx+3*h/4,ww)
   ww = yy - (3*k1 - 2*k2 - 12*k3 + 12*k4 - 8*k5) / 7
   k6 = h * f(xx+h,ww)
   yy = yy + (7*k1 + 32*k3 + 12*k4 + 32*k5 + 7*k6) / 90
   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 rk5.dat")