# rk1: one step in the standard: Runge-Kutta  method for the ODE
#
#   dy/dx = f(x,y)
#
# with initial conditions y(x0) = y0.   (see also file "desolvers")
#
# Example:
#
#   f := (x,y) -> y
#   rk1( 0.0, 1.0, 0.1, f);
#
# Usage:
#
#   rk( x0, y0, n, f );
#

rk1:=proc(x0,y0,h,f)
local k1,k2,k3,k4,s:
 k1:=h*f(x0,y0):
 k2:=h*f((x0+h/2),(y0+k1/2)):
 k3:=h*f((x0+h/2),(y0+k2/2)):
 k4:=h*f((x0+h),(y0+k3)):
 s:=y0+(k1+2*k2+2*k3+k4)/6:
 RETURN(s):
end:


# rk: repetitively executes rk1
#
# Example:
#
#   f := (x,y) -> y
#   rk( 0.0, 1.0, 1.0, 100, f )
#
# This computes an approximation to $e$ by solving
# dy/dx = y with y(0.0) = 1.0 using the Runge-Kutta
# method with 100 subdvisions.  Note that to
# compute $e$ using Maple, do this:
#
#   evalf(E);
#
# Usage:
#
#  rk( a, b, Y, n, f)
#
# The value computed is an approximation to y(b) where
# y(x) is a solution to dy/ds = f(x,y) with y(a) = Y.
#


rk := proc( x0, x1, y0, n, f)

  local x, y, h, i:               # set up internal variables
  h := evalf((x1 - x0)/n):        # compute step size
  x := x0; y := y0:               # initialize

  for i from 1 to n do
    y := rk1( x, y, h, f ):       # compute new y-value
    x := x + h:                   # compute new x-value
  od:

  RETURN(y);                      # value of rk(...)
end:


rkplot := proc( x0, x1, y0, n, f)
      local i, x, y, h, data:         # set up internal variables
      h := evalf((x1 - x0)/n):        # compute step size
      x := x0; y := y0:               # initial point  
      data := array(0..n);            # array to hold the data
      data[0] := [x0,y0];
      for i from 1 to n do
        y := rk1( x, y, h, f ):       # compute new y-value
        x := x + h:                   # compute new x-value
        data[i] := [x,y]
      od:  
      plot( convert( data, list ) );  # convert array to list & plot
    end: