# Heun1: one step in Heun's method for the ODE
#
#   dy/dx = f(x,y)
#
# with initial conditions y(x0) = y0.   (see also file "desolvers")
#
# Example:
#
#   f := (x,y) -> y
#   Heun1( 0.0, 1.0, 0.1, f);
#
# Usage:
#
#   Heun( x0, y0, h, f );
#

Heun1:=proc(x0,y0,h,f)
  local s,s0:
   s0 := f(x0,y0):
   s := y0+h*(s0+f((x0+h),(y0+h*s0)))/2:
 RETURN(s):
end:


# Heun: repetitively executes Heun1
#
# Example:
#
#   f := (x,y) -> y
#   Heun( 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 Heun's method
# with 100 subintervals.  Note that to
# compute $e$ using Maple, do this:
#
#   evalf(E);
#
# Usage:
#
#  Heun( a, b, Y, h, 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.
#


Heun := 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 := Heun1( x, y, h, f ):     # compute new y-value
    x := x + h:                   # compute new x-value
  od:

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


Heunplot := 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 := Heun1( 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: