# Numerical Solution of ODE - Euler's Method and Improved Euler 
#             ODE- IVP :   y' = f(x,y),   y(0) = 1 
#  
# Goal:  Construct graphs comparing the exact solution to the numerical 
#          approximations. 
#  
# Turn in graphs for the following two IVPs that illustrate the use of 
# two numerical methods and the exact solution. 
#  
#         a.   f(x,y) = y*(2-y) 
#  
#         b.   f(x,y) = x- y^2
> restart;
# The sequence of approximations is (x[n],y[n]),  n=0,1,2,... 
#  
# Define initial conditions and step size 
> x[0]:= 0; y[0] := 1.0; h:= 0.1;
# Define the function f(x,y) in the RHS of ODE 
> f := (x,y) -> y;
# Euler's method or algorithm:
>
for n from 1 to 20 do\
  x[n] := n*h;\
  y[n] := y[n-1] + h*f(x[n-1],y[n-1]);\
od:
# The next command generates sequence of of approximatons:
> data := [seq([x[n],y[n]],n=0..20)]:
# Generate plot which is not displayed but instead stored under the name t1
> t1 := plot(data,style=point,color=red):
# Improved Euler's method -  use different names for approximation
> xx[0]:=0: yy[0]:=1:\
for n from 1 to 20 do\
  xx[n] := n*h;\
  ystar := y[n-1] + h*f(x[n-1],y[n-1]);\
  yy[n] := yy[n-1] + \
  h/2.0*(f(xx[n-1],yy[n-1])+f(xx[n],ystar)):\
od:
# Generate the sequence of approximations for the Improved Euler method
> data_improve := [seq([xx[n],yy[n]],n=0..20)]:
# Generate plot which is not displayed but instead stored under the name t3
> t3 := plot(data_improve,style=point,color=blue):
# Now have Maple construct the exact solution, if possible. 
# We write symbolic form of ODE using a new function u(x)
> eqn := diff(u(z),z) = f(z,u(z));
# Try to construct exact solution of IVP using dsolve
> dsolve({eqn,u(0)=1},u(z));
# We need to define u as a function
> u := unapply(rhs(%),z);
# Graph without displaying the exact solution
> t2 := plot(u(z),z=0..2):
# Plot three graphs on same axes but insert your name as the title.
> plots[display]({t1,t2,t3},title=`your name`);
>