> 
# Euler's Method
# The Picard proof-idea  has the drawback of requiring a large amount of
# symbolic integration. What is worse, in many cases it is not possible
# to work the integrals that arise. So it has more theoretic than
# numerical use. Here we give another perspective on the business of
# "why are there solutions"?
# 
# Euler's method, like the Picard idea, has to do with solutions to
# first order differential equations. Again, we will illustrate it on
# the differential equation y'=x+y^2. 
> f:=(x,y)->x+y^2;

                                            2
                        f := (x, y) -> x + y

> x0:=0;y0:=0;

                               x0 := 0


                               y0 := 0

> pointsequence:=NULL;

                          pointsequence :=

> dx:=0.1;

                               dx := .1

> x:=x0;y:=y0;for k from 1 to 20 do
> dy:=dx*f(x,y):x:=x+dx:y:=y+dy:pointsequence:=pointsequence,[x,y]
> od:pointlist:=[pointsequence]:plot(pointlist);

                                x := 0


                                y := 0


> pl1:=":
> 
> pointsequence:=NULL;dx:=0.05;x:=x0;y:=y0;for k from 1 to 40 do
> dy:=dx*f(x,y):x:=x+dx:y:=y+dy:pointsequence:=pointsequence,[x,y]
> od:pointlist:=[pointsequence]:plot(pointlist);
> 

                          pointsequence :=


                              dx := .05


                                x := 0


                                y := 0


> pl2:=":
> with(plots);

[animate, animate3d, changecoords, complexplot, complexplot3d,

    conformal, contourplot, contourplot3d, coordplot, coordplot3d,

    cylinderplot, densityplot, display, display3d, fieldplot,

    fieldplot3d, gradplot, gradplot3d, implicitplot, implicitplot3d,

    inequal, listcontplot, listcontplot3d, listdensityplot, listplot,

    listplot3d, loglogplot, logplot, matrixplot, odeplot, pareto,

    pointplot, pointplot3d, polarplot, polygonplot, polygonplot3d,

    polyhedraplot, replot, rootlocus, semilogplot, setoptions,

    setoptions3d, spacecurve, sparsematrixplot, sphereplot, surfdata,

    textplot, textplot3d, tubeplot]

> display(pl1,pl2);

> 
# Well, this is certainly a simpler idea. You just follow the direction
# field out along a straight line of slope m=f(x_k,y_k) for a horizontal
# distance dx, changing x to x+dx and y to y+m*dx=y+dy. Then connect the
# dots. 
# 
# Is it a perfect idea? It has half of what computer scientist David
# Gelernter calls machine beauty. It has simplicity. But machine beauty
# is a blend of simplicity and power, and as numerical methods go, this
# one is poor. Too much arithmetic must be done for too little accuracy.
# 
# 
# There are all kinds of things that can be done to improve on this
# basic idea, and there are some elegant algorithms. Rather than delve
# into the details, suffice it to say that these refinements of Euler's
# idea have been incorporated into Maple. Thus, we can ask it for a
# numerical solution to the differential equation and if we don't botch
# the syntax, we'll get an answer. 
> with(DEtools);

[DEnormal, DEplot, DEplot3d, Dchangevar, PDEchangecoords, PDEplot,

    autonomous, convertAlg, convertsys, dfieldplot, indicialeq,

    phaseportrait, reduceOrder, regularsp, translate, untranslate,

    varparam]

> dsolve({diff(y(x),x)=x+y^2,y(0)=0},type=numeric);
Error, (in dsolve) dsolve/inputck1 expects its 2nd argument, vars, to
be of type {list, set, function}, but received type = numeric
> dsolve({diff(y(x),x)=x+y^2,y(0)=0},y(x),type=numeric);

                       proc(rkf45_x)  ...  end

> whazzat:=";

                  whazzat := proc(rkf45_x)  ...  end

> whazzat(1);

                  [x = 1, y(x) = .5571617684443275]

# This is unfortunate, but Maple's numerical solutions come to us in a
# format that is inconvenient for the use we have in mind next. We shall
# have to extract the right hand side of the second part of the list
# above, to get at the number 0.557161... and we shall have to generate
# a list of points so we can plot them. 
# 
# (Maple has another command, DEplot, which shortcuts all that. But
# right now we're looking at the numerical dsolver.)
> ptlist3:=[seq([k/10,rhs(whazzat(k/10)[2])],k=0..19)];

ptlist3 := [[0, 0], [1/10, .005000500074847118],

    [1/5, .02001601607101244], [3/10, .04512191161064932],

    [2/5, .08051612991663130], [1/2, .1265873093237576],

    [3/5, .1839959448121971], [7/10, .2537802508902369],

    [4/5, .3375056819032614], [9/10, .4374901365628385],

                             11
    [1, .5571617579371073], [--, .7016565745079046],
                             10

                               13
    [6/5, .8788712732258639], [--, 1.101436407988666],
                               10

    [7/5, 1.390709325772345], [3/2, 1.785693450133243],

                               17
    [8/5, 2.365808454820888], [--, 3.321367354223440],
                               10

                               19
    [9/5, 5.250957094058988], [--, 11.52504961002434]]
                               10

> pl3:=plot(ptlist3):
> display(pl1,pl2,pl3);

# So the appeal to the built-in Maple numerical solver shows a rather
# striking divergence from the Euler's method with step sizes 1/10 or
# 1/20, even by x=1.5. We would have to take a very small step size, and
# pay by taking very many steps, if we wanted to get comparable accuracy
# out of Euler's method.