User must enter: function f(x,y), Nsteps, a, b, y0 below.
This code is set up to use equally spaced nodes (constant stepsize h).

> restart ;    (This clears Maple's memory. Always best to make this the first line in a Maple code.)

> Digits := 16 ;    (Sets the precision of Maple. This should be an integer.)

> f := ( x, y ) –> evalf( y * cos(x) ) ;
    (This defines f as a function of x and y. You may use any variable names you want, but it's best not to use names you've used elsewhere.)

> a := 0.0 ;   b := 2.0 ;    (We are solving the problem on the interval [0.0, 2.0]. a is x₀ from the initial condition.)

> y0 := 1.0 ;

> Nsteps := 10 ;    (We are using 10 steps. This should be an integer.)

Euler's method follows:

> h := (b - a) / Nsteps ;    (Determine (set) the nodal stepsize.)

> X[0] := a ;   Y[0] := y0 ;
    ( Use square brackets, not parentheses. Square brackets are used to denote elements of a vector or elements of a sequence. )

> printf("\n                x                   y\n") :    (after each line in this block, hit Shift-Enter)
    printf("-----------------------------------------------\n") :
    printf("  %3d  %18.12f  %18.12f\n", 0, X[0], Y[0] ) :    (This prints the initial value.)
    for n from 0 to (Nsteps–1) do
          X[n+1]   :=   a + (n+1) * h :
          Y[n+1]   :=   Y[n] + h * f( X[n], Y[n] ) :
          printf("  %3d  %18.12f  %18.12f\n", n+1, X[n+1], Y[n+1] ) :
    od:

   n             x                   y
-----------------------------------------------
    0      0.000000000000      1.000000000000
    1      0.200000000000      1.200000000000
    2      0.400000000000      1.435215978682
    3      0.600000000000      1.699600269869
    4      0.800000000000      1.980148396635
    5      1.000000000000      2.256064931323
    6      1.200000000000      2.499856348239
    7      1.400000000000      2.681024814812
    8      1.600000000000      2.772162040375
    9      1.800000000000      2.755972878911
   10      2.000000000000      2.630740316710

> plot( [ [ X[k], Y[k] ]$k = 0 .. Nsteps ], a .. b, 0.0 .. 3.0, style = point, symbolsize = 30, symbol = diamond ) ;

> plot( [ [ X[k], Y[k] ]$k = 0 .. Nsteps ], a .. b, 0.0 .. 3.0, style = line, thickness = 4 ) ;