#
# Solution to the problem:
#   y' = f(x,y)
# by using f(x,y) to determine the direction of the solution from an 
# initial value.  That is, a graphical represention of Euler's Method.
# 
#                              Maple Midterm Project #1
#                                      Euler's Method
# 
# (Solving a differential equation via the direction field.)
# 
--------------------------------------------------------------------------------
# 
# First, re-initialize all variables in Maple...
# 
> restart;\

--------------------------------------------------------------------------------
# 
# Part a) Write a Maple Program which takes a given function f, 
# initial conditions, and h,
# and then generates a sequence of points, 
# which it then plots.
# 
--------------------------------------------------------------------------------
# 
# Given a differential equation:
# 
> diff ( y(x), x) = f(x,y);

                                 d
                               ---- y(x) = f(x, y)
                                dx
# and initial conditions:
# 
> y(x0) = y0;

                                   y(x0) = y0
--------------------------------------------------------------------------------
# 
# Define the right hand side of the equation:
# 
> f := proc (x,y)\
          y;\
      end:\

--------------------------------------------------------------------------------
# 
# Define the number of points, N, 
# 
> N := 100:\

# and the stepsize, h:
# 
> h := 0.01:\

--------------------------------------------------------------------------------
# 
# Define the initial conditions
# x0 and y0:
# 
> x0 := 0:\
y0 := 1:\

--------------------------------------------------------------------------------
# 
# Set x[0] and y[0] to be the initial conditions:
# 
> x[0] := x0:\
y[0] := y0:\

--------------------------------------------------------------------------------
# 
# Now begin loop for 
# x[n] and y[n]:
# 
> for n from 0 to N do\
   x[n+1] := x[n] + h;\
   y[n+1] := y[n] + h*f(x[n],y[n]);\
od:\

--------------------------------------------------------------------------------
# 
# Finally, plot the points:
# 
> pts := [ [ x[i], y[i] ] $i=0..N ]:\
plot(pts);\

--------------------------------------------------------------------------------
# 
# Part b) Calculate the error, and investigate how it depends on initial
# conditions, f, and h.
# 
# Define the solution, yexact,
# 
> yexact := x -> exp(x);

                                  yexact := exp
--------------------------------------------------------------------------------
# 
# Calculate the error:
# 
> for n from 0 to N do\
   error[n] := yexact(x[n]) - y[n];\
od:
# 
--------------------------------------------------------------------------------
# 
# Graph the error, vs. x:
# 
> err := [ [ x[i], error[i] ] $i=0..N ]:\
plot(err,title=`error vs x`);\

--------------------------------------------------------------------------------
# 
# Graph the log of the error, vs. x:
# 
> logerr := [ [ x[i], log(abs(error[i])) ] $i=1..N ]:\
plot(logerr,title=`log error vs x`);\

--------------------------------------------------------------------------------