-----------------------------------------------------
Example of Maple Code 	(by Vladimir Dobrushkin)
-----------------------------------------------------

> restart;
> with(DEtools):
> with(plots):

## define a function: 
> f:=(x,y) -> x^2 - y^2; ## for example;

## and the initial condition and step size (h): 
> x[0] := x0: y[0] := y0: h:=0.01: 

## define the terminated point:
> xn := xn;
## define the number of points: 
> n := (xn-x0)/h;
> nm :=n-1;

#### A) Euler method:

> for k from 0 to nm do
	k1 := f(x[k], y[k]);
	y[k+1] := y[k] + h*k1;
	x[k+1] := x[k] + h;
> od:
## The Euler method can be implemented using subrouting "proc."
## Let us consider a Cauchy problem for a differential equation 
## y' = f(x,y),  y(a) =A   ## on an interval [a,b]

## Let N be the number of steps required to reach the point b

> eu:=proc(f,a,b,A,N)
  local n,xn, yn, h, k1, w;
      h:=evalf((b-a)/N); xn:=evalf(a); yn:=evalf(A); 
      w[0]:=[xn,yn];
          for n from 1 to N do 
          k1:=evalf(f(xn,yn));
          yn:=evalf(yn+h*k1); 
          xn:=evalf(xn+h); 
          w[n]:=[xn,yn];
          end do;
      w; 
   end;

## To call the procedure  for a given say a=0, b=1, A=0, and N=10
   type: 

>  u:=eu(f,0,1,0,10);

## To see the yn (n-th value of y) at point, say 6, type: 

> u[6][2]; 

## To print values xn (=a+n*h) and yn type: 

> for n from 0 to N do print(u[n][1], u[n][2]) od; 
## or 
> for n from 0 to N do print(x||n=u[n][1],y||n=u[n][2]) od; 

## If you want to decrease the stepsize h in 10 times: 

> u2:=eu((f,0,1,0,N*10);
> for n from 0 to N*10 by 10 do 
  print(x||n=u[n][1],y||n=u[n][2]) od; 


## define the sequence of points: 
> pts:=seq([x[k],y[k]],k=0..n);
# or using proc 
> pts1:=seq(u[k],k=0..N);

## and plot them: 
> pointplot({seq([x[k],y[k]],k=1..n)});


#### B) Runge-Kutta method: 

> for k from 0 to nm do   
    k1:=f(x[k],y[k]);
    k2:=f(x[k]+h/2,y[k]+k1*h/2);
    k3:=f(x[k]+h,y[k]+2*h*k2-h*k1);
    x[k+1]:=x[k]+h;
    y[k+1]:=y[k]+(h/6)*(k1+4*k2+k3);
 od:
> pts:=seq([x[k],y[k]],k=0..n);





#### Exact solution (for ODE y' = x^2 + y^2):
> dsolve({diff(y(x),x)=y(x)^2+x^2,y(0)=0},y(x));
> fsol:=x -> -x*(-BesselJ(-3/4,1/2*x^2)+BesselY(-3/4,1/2*x^2))/
              (-BesselJ(1/4,1/2*x^2)+BesselY(1/4,1/2*x^2)):
> psol:=plot(fsol(x),x=0..1.5):
> display(psol);

#Example 1.b ------------------------------------------------
> x[0]:=0.:  y[0]:=0: h := 0.15:
> for k from 0 to 9 do
    x[k+1] := x[k]+h:
    y[k+1] := y[k]+(x[k]^2+y[k]^2)*h:
  od:
> ptsapp1 := seq([x[k],y[k]],k=1..10):
> ptssol1 := seq([x[k],fsol(x[k])],k=1..10):
> p1app:=pointplot({ptsapp1}):
> display(psol,p1app);
> errorplot1 := pointplot( {seq([x[k],fsol(x[k])-y[k]],k=1..10)}):
> display(errorplot1);

#Example 2.a:-------------------------------------------------
> x[0]:=0.:  y[0]:=0: h := 0.015:
> for k from 0 to 99 do
    x[k+1] := x[k]+h:
    y[k+1] := y[k]+(x[k]^2+y[k]^2)*h:
  od:
> ptsapp2 := seq([x[k],y[k]],k=1..100):
> ptssol2 := seq([x[k],fsol(x[k])],k=1..100):
> p2app:=pointplot({ptsapp2}):
> display(psol,p2app);
> errorplot2 := pointplot( {seq([x[k],fsol(x[k])-y[k]],k=1..100)}):
> display(errorplot2);

#Example 2.b -------------------------------------------------
> x[0]:=0.:  y[0]:=0: h := 0.0015:
> for k from 0 to 999 do
    x[k+1] := x[k]+h:
    y[k+1] := y[k]+(x[k]^2+y[k]^2)*h:
  od:
> ptsapp3 := seq([x[k],y[k]],k=1..1000):
> ptssol3 := seq([x[k],fsol(x[k])],k=1..1000):
> p3app:=pointplot({ptsapp3}):
> display(psol,p3app);
> errorplot3 := pointplot( {seq([x[k],fsol(x[k])-y[k]],k=1..1000)}):
> display(errorplot3);

# Example 2.c --------------------------------------------------
> display(errorplot1,errorplot2,errorplot3);
> errorplot12c := logplot( [seq([x[k],fsol(x[k])-y[k]],k=1..100)],symbol=box):
> display(errorplot12c);

# Example 2 ------------------------------------------------------
> x[0]:=0.:  y[0]:=0: h := 0.15:
> for k from 0 to 9 do
    x[k+1] := x[k]+h:
    y[k+1] := y[k]+(x[k]^2+y[k]^2)*h + (2*x[k] + 2*y[k]*(x[k]^2+y[k]^2) )*h^2/2:
  od:
> ptsapp4 := seq([x[k],y[k]],k=1..10):
> ptssol4 := seq([x[k],fsol(x[k])],k=1..10):
> p4app:=pointplot({ptsapp4}):
> display(psol,p4app);
> errorplot4 := pointplot( {seq([x[k],fsol(x[k])-y[k]],k=1..10)}):
> display(errorplot4);  #<additional commands not showing because of lack of space>