@sarra 

Thank you, error is a maple word. It's changed.

I have a question, we solve the Ode, in the interval [0,1] but the aprroximate solution is given on some points 0,..,0.22 in your code send me this morning.

Please see the hole code:

First problem: Why the discretisation take not all the point of the interval, ideed we use x[n+1]=x[n]+h;

                                    NULL
restart:
 
ode := diff(y(x), x) = -y(x);
f:=(x,y)->-y;
analyticsol := rhs(dsolve({ode, y(0) = 1}));
## Procedure adaptative RK23
RKadaptivestepsize := proc (f, a, b, epsilon, N)
local n, x, y,k,z,R,p,h,hstar,prec1;
if nargs<2 then
      error "invalid arguments"
   end if;
p:=2;
h := evalf(b-a)/N;
x[0] := a; y[0] := 1;z[0]:=1;
for n from 0 to N-1 do
x[n+1] := a+(n+1)*h;

k[1] := f(x[n], y[n]);
k[2] := f(x[n]+h, y[n]+h*k[1]);
k[3] := f(x[n]+h/2, y[n]+h/4*(k[1]+k[2]));
## 2-stage runge Kutta.
z[n+1] := z[n]+(h/2)*(k[1]+k[2]);
y[n+1] := y[n]+(h/6)*(k[1]+k[2]+4*k[3]);

## local Error
R:=abs(y[n+1]-z[n+1]);
hstar:=sqrt(epsilon/R);
if R<=epsilon then
 h:=h; # we take the same stepsize
 x[n+1] := x[n]+h;
 y[n]:=y[n+1];
 z[n]:=z[n+1];
 n:=n+1;
else
h:=hstar;
end if;
end do;
return [seq([x[i], z[i]], i = 0 .. N)];
#[seq([x[i], z[i]], i = 0 .. N)];
end proc:
epsilon:=1e-7;

### compute the error
RKadaptivestepsize(f,0,1,epsilon,20);  ## we work in [0,1] but the approximate solution is given from 0, to  0.5471001066, see please the next table.


                                   
[[0, 1], [0.05000000000, 0.9512500000], [0.1385641755, 0.8876298110],

     [0.1106146866, 0.8555075471], [0.1405708399, 0.8259737692],

     [0.1693613798, 0.7984728927], [0.1972701478, 0.7726549758],

     [0.2244820369, 0.7482771818], [0.2511259696, 0.7251600066],

     [0.2772963821, 0.7031647053], [0.3030646746, 0.6821804387],

     [0.3284870239, 0.6621163263], [0.3536081857, 0.6428963231],

     [0.3784648380, 0.6244557160], [0.4030871298, 0.6067386783],

     [0.4275008760, 0.5896964621], [0.4517273038, 0.5732861030],

     [0.4757852295, 0.5574693832], [0.4996908733, 0.5422120558],

     [0.5234582872, 0.5274832335], [0.5471001066, 0.5132548894]]

## here we compute the error at x=1, between the numerical obtained in previous ( must be corrcted) ##and the exact solution.

dataerror:= N->evalf(abs(RKadaptivestepsize(f,0,1,epsilon,N)[1+N][2]-(eval(analyticsol, x = 1))));

## Error between the true solution and Numerical solution obtainted from ##RKadaptivestepsize
for N from 2 by 2 to 500 do
dataerror:= N->evalf(abs(RKadaptivestepsize(f,0,1,epsilon,N)[1+N][2]-(eval(analyticsol, x = 1))));
##  sequence of data error
data[err] := [seq([N, dataerror(N)], N = 2 .. 500, 2)]:
if  data[err][N][2]<=epsilon then   
printf("%a  is the number of steps required using 3-step Runge Kutta Method to achieve an  eroor of 1e-6 .", N)
end if;   
end do;
end do;







        Thank your for your help. (merci d'avance pour vos remarques).

@sarra 

This link contain some methods.  thank for helping me.

http://www.peterstone.name/Maplepgs/numDE.html

@sarra 

I find this code I whish That my code is correct. What do you think?
comparewithfcn := proc(pts,f)
   local xx,y1,y2,r,i,j,n,rows,fn,x,prec,prcsn,saveDigits,
   startoptions,Options,md,t,maxerr,format1,format2,g,e,zero,
   proctype,ertyp,xmax,prec1,prec2,vars,averr;

   proctype := false;
   if nargs<2 then
      error "invalid arguments; the basic syntax is 'comparewithfcn([[x1,y1],..,[xn,yn]],f(x),x)' or 'comparewithfcn([[x1,y1],..,[xn,yn]],f)'"
   end if;
   if not type(pts,listlist) then
      error "the 1st argument, %1, is invalid .. it should be a list of points",pts;
   end if;
   if nargs>2 and not type(args[3],equation) then
      x := args[3];
      if not type(x,name) then
         error "the 3rd optional argument must be the name of the independent variable"
      end if;
      startoptions := 4;
      vars := indets(f,name) minus indets(f,realcons);
      if vars<>{x} then
         if not type(f,algebraic) or not has(indets(f),{Int,Sum}) then
            error "the 2nd argument, %1, is invalid .. it should be an expression which depends only on the single variable %2",f,x;
         end if;
      end if;
   else
      if type(f,procedure) or (type(f,`@`) and type({op(f)},set(procedure))) then
         proctype := true;
         startoptions := 3;
      else
         error "the 2nd argument, %1, is invalid .. it should be a function of one variable or an expression in the one variable given as a 3rd argument",f;
      end if;
   end if;   

   # Get the options.
   md := 0;
   ertyp := 1;
   if nargs>=startoptions then
      Options :=[args[startoptions..nargs]];
      if not type(Options,list(equation)) then
         error "each optional argument must be an equation"
      end if;
      if hasoption(Options,'mode','md','Options') then
         if not (md='linebyline' or md='matrix') then
            error "\"mode\" must be 'matrix' or 'linebyline'"
         end if;
         if md='linebyline' then md := 0 else md := 1 end if;
      end if;
      if hasoption(Options,'errtype','ertyp','Options') or
         hasoption(Options,'errortype','ertyp','Options') then
         if not member(ertyp,{'absolute','ABSOLUTE','ABS','relative','RELATIVE','REL'}) then
            error "\"errtype\" option must be 'absolute' <-> 'ABSOLUTE' <-> 'ABS' or 'relative' <-> 'RELATIVE' <-> 'REL'"
         end if;
         if member(ertyp,{'absolute','ABSOLUTE','ABS'})
         then ertyp := 0 else ertyp := 1 end if;
      end if;
      if nops(Options)>0 then
         error "%1 is not a valid option for %2",op(1,Options), procname;
      end if;
   end if;
   n := nops(pts);
   prcsn := 0;

   # Check the data and find its maximum precision.
   for i to n do
      if nops(pts[i])<>2 then
         error "the 1st argument must be a list of points, where each point is itself a list with two members"
      end if;
      t := pts[i,1];
      if type(t,float) and type(t,numeric) then
         prec1 := length(convert(op(1,t),string));
      elif type(t,realcons) then
         prec1 := Digits;
      else
         error "the 1st argument must be a list of points, where each point is itself a list of two real numbers"
      end if;
      t := pts[i,2];
      if type(t,float) and type(t,numeric) then
         prec2 := length(op(1,t));
      elif type(t,realcons) then
         prec2 := Digits;
      else
         error "the 1st argument must be a list of points, where each point is itself a list of two real numbers"
      end if;
      prec := max(prec1,prec2);
      if prec>prcsn then prcsn := prec end if;
   end do;
   saveDigits := Digits;
   Digits := prcsn;
   prec := trunc(prcsn/2);
   
   if proctype then
      fn := f;
   else
      fn := unapply(evalf(f),x);
   end if;
   rows := NULL;
   maxerr := 0;
   averr := 0;
   zero := false;
   format1 := cat("%",convert(prcsn+5,string),".",convert(prcsn,string),g);
   format2 := cat("%",convert(prec+4,string),".",convert(prec-1,string),e);

   for i from 1 to n do
      xx := evalf(pts[i,1]);
      y1 := evalf(pts[i,2]);
      y2 := traperror(evalf(fn(xx)));
      if y2=lasterror or not type(y2,numeric) then
         error "function failed to evaluate to a real floating point number at %1",xx;
      end if;
      
      if ertyp=0 then
         r := evalf(abs(y1-y2));
         averr := averr+r;
         if r>maxerr then
            maxerr := r;
            xmax := xx;
         end if;
      else
         if y2<>0 then
            r := evalf(abs(y1-y2)/abs(y2));
            if r>maxerr then
               maxerr := r;
               xmax := xx;
            end if;
         else
            zero := true;
            r := infinity;
         end if;
      end if;

      if md=0 then
         printf(format1,xx);
         printf(` `);
         printf(format1,y1);
         printf(`   function val: `);
         printf(format1,y2);
         if ertyp=0 then
            printf(`     abs err: `);
            printf(format2,r);
         else
            printf(`     rel err: `);
            if y2<>0 then
               printf(format2,r)
            else
               printf(infinity);
            end if;
         end if;
         printf(`\n`);
      else
         rows := rows,[xx,y1,y2,r];
      end if;
   end do;
   
   print(``);
   if ertyp=0 then
      printf(`              Maximum absolute error: `);
   else
      printf(`              Maximum relative error: `);      
   end if;
   printf(format2,maxerr);

   if maxerr<>0 then
      printf(`\n              obtained for the input value: `);
      printf(format1,xmax);
   end if;
   
   if ertyp=1 and zero then
      printf(`\n              excluding any cases where the function value is zero.`);
   end if;

   if ertyp=0 then
      averr := averr/n;
      print(``);
      printf(`              Mean absolute error: `);
      printf(format2,averr);    
   end if;

   Digits := saveDigits;
   if md=1 then
      print(``);
      if ertyp=0 then
         return array([[x,"discrete value","function value","absolute err"],rows]);
      else
         return array([[x,"discrete value","function value","relative err"],rows]);
      end if;
   else
      return NULL;
   end if;
end proc:

@mehdi jafari 

Thank your for tour astutute remarks.

According to the previous code. Using it, one can obtained an approximate solution to the ode.

Now, I would like to compute the error between the approximate solution abtained and the axact solution.
I add these lines, I feel that they are correct but there is an error when I make RUN. 

In order to compute the error between the Numerical solution and exact. Then I say when the error is of order epsilon break and display the "n" ( the step n which give the error epsilon).


restart:
ode := diff(y(x), x) = -y(x);
f:=(x,y)->-y;
analyticsol := rhs(dsolve({ode, y(0) = 1}));
## Procedure adaptative RK23

RKadaptivestepsize := proc (f, a, b, epsilon, N)
local n, x, y,k,z,R,p,h,hstar;
p:=2;
h := evalf(b-a)/N;
x[0] := a; y[0] := 1;z[0]:=1;
for n from 0 to N-1 do
x[n+1] := a+(n+1)*h;
k[1] := f(x[n], y[n]);
k[2] := f(x[n]+h, y[n]+h*k[1]);
k[3] := f(x[n]+h/2, y[n]+h/4*(k[1]+k[2]));
## 2-stage runge Kutta.
z[n+1] := z[n]+(h/2)*(k[1]+k[2]);
y[n+1] := y[n]+(h/6)*(k[1]+k[2]+4*k[3]);

## local Error
R:=abs(y[n+1]-z[n+1]);
hstar:=sqrt(epsilon/R);
if R<=epsilon then
 h:=h; # we take the same stepsize
 x[n+1] := x[n]+h;
 y[n]:=y[n+1];
 z[n]:=z[n+1];
 n:=n+1;
else
h:=hstar;
end if;
end do;
return [seq([x[i], z[i]], i = 0 .. N)];
#[seq([x[i], z[i]], i = 0 .. N)];
end proc:
epsilon:=1e-6;
dataerror:= N->evalf(abs(RKadaptivestepsize(f,0,1,epsilon,N)[1+N][2]-(eval(analyticsol, x = 1))));
############### New lines added.
## Error between the true solution and Numerical solution obtainted from ##RKadaptivestepsize
## here, I compute the error
for N from 2 by 2 to 500 do
Test_error:= N->evalf(abs(RKadaptivestepsize(f,0,1,epsilon,N)[1+N][2]-(eval(analyticsol, x = 1))));
##  sequence of data error
data[error] := [seq([N, Test_error(N)], N = 2 .. 500, 2)]:
if  data[error][k][2]<=epsilon then   
printf("%a  is the number of steps required using 3-step Runge Kutta Method to achieve an  eroor of 1e-6 .", k)
break ;
end if;   
end do;
end do;







                 
                          
Error, reserved word `error` unexpected

@sarra 

How can add this sentence and get an answer from my code. one again thank you for your remarks

printf("%a  is the number of steps required using 3-step Runge Kutta Method to achieve an  eroor of 10^(-6) .", n);

@mehdi jafari 

Thank you

Thank you for your remark
Here we use z numerical solution 2-stage Runge Kutta, and y numerical solution  of 3-stage runge Kutta.

I try this method but I have a problem to run the code
 

ode := diff(y(x), x) = x-y(x);
f:=(x,y)->-y;
analyticsol := rhs(dsolve({ode, y(0) = 1}));

## Numerical solution
epsilon:=1e-5;
RKadaptivestepsize := proc (f, a, b, epsilon, N)
local x, y, n, h,k,z,R,p;
p:=2;
h := evalf(b-a)/N; ## we begin with this setpsize
x[0] := a; y[0] := 1; ## Initialisation
for n from 0 to N-1 do  ##loop
x[n+1] := a+(n+1)*h;  ## noeuds
k[1] := f(x[n], y[n]);
k[2] := f(x[n]+h, y[n]+h*k[1]);
k[3] := f(x[n]+h/2, y[n]+h/4*(k[1]+k[2]));
z[n+1] := z[n]+(h/2)*(k[1]+k[2]);## 2-stage runge Kutta.
y[n+1] := y[n]+(h/6)*(k[1]+k[2]+4*k[3]);
R:=abs(y[n+1]-z[n+1]);
hstar:=h*sqrt(epsilon/R);
if R<epsilon    then
     h:=1.5*h;
    #x[n] := x[n+1]+h;
     return z[n];
print(n,x[n])
else
h:=0.4*hstar;
n:=n+1;
end if
 end do;
[seq([x[n], y[n]], n = 0 .. N)];
[seq([x[n], z[n]], n = 0 .. N)];
end proc:

@Adri van der Meer 

Thank you.

Are you agree about the definition of N using round...

N:=round(x0/(2*h));

for me it's seem to get some point outside the interval of computation, [0,x0].

Mybe must I use floor instead of round.

@Carl Love 

I agree with you.

Here, I my case i can expect that :

Forward Euler method :  epsilon(h):=O(h^2);

Huen Method :  epsilon(h):=O(h^3));

4-stage Range Kutta   epsilon(h):=O(h^5);

what do you think about these results.

@Carl Love 

Thank you.  I try yours, it's working. I think there is a command maple using isolate... but how.... 

In the case where A= [[a,a],[2*a,2*a]];

how can we get using maple A:=a*[[1,1],[1,2]];

@Joe Riel 

Thanks you for your remarks

@sarra 

I want to plot the error defined in the last line.  I tried my code, Heun1 work well, I dont konw why Heun2 don't working. Please, I need a Help.

here the code with maple code, without pciture.

f :=(x,y)-> 1/(1+sin(y));

ode:=diff(y(x),x)=f(x,y);

## Heun Method with stepsize 2*h

Heun1 := proc(f, x0, h)
local x, y, i, N, k;
N := round((1/2)*x0/h);
y := Array(0 .. N);
x := Array(0 .. N);
x[0] := 0;
y[0] := (1/4)*Pi;
for i from 0 to N-1 do
x[i+1] := (i+1)*2*h;
k[1] := f(x[i], y[i]);
k[2] := f(x[i]+h, y[i]+h*k[1]);
y[i+1] := y[i]+(1/2)*h*(k[1]+k[2]);
 end do;
 evalf([seq([x[i], y[i]], i = 0..N)])
 end proc

## Heun Method with stepsize h

Heun2 := proc(f, x0, h)
local x, y, i, N, k;
N := round((1/2)*x0/h);
y := Array(0 .. 2*N);
x := Array(0 .. 2*N);
x[0] := 0;
y[0] := (1/4)*Pi;
for i from 0 to 2*N-1 do
x[i+1] := (i+1)*h;
k[1] := f(x[i], y[i]);
k[2] := f(x[i]+h, y[i]+h*k[1]);
y[i+1] := y[i]+(1/2)*h*(k[1]+k[2]);
 end do;
 evalf([seq([x[2*i], y[2*i]], i = 0..N)])
 end proc

 Heun1((x,y)->f(x, y), 1,1e-2);  ## Heun1 work very well

 Heun2((x,y)->f(x, y), 1,1e-2); ## Dont working...

## The error

x0:=1;

epsilon:= (x0,h)->sqrt(1/(( round(x0/(2*h))+1))*add(( Heun1((x,y)->f(x,y),x0,round(x0/(2*h)))[i][2]-Heun2((x,y)->f(x,y),x0,  round(x0/(2*h)))[i][2])^(2),i=1..N+1));

data[Heun] := [seq([h,epsilon(x0,h)], h =1e-2..1)];

loglog(epsilon(1,h),h=1e-2..1);

@mehdi jafari 

Sorry to put the code as a picture. The code was modified according to your remark.

I hope get a solution. Thank you.

@Joe Riel 

Thank you for your help.  It is okay.

@Joe Riel 

Thank you for your answer. 

Mu code RK3 (f,a,b, y0,N) work well. 

Here I would libe to make a loop on "N" number of step to get error about 0.001. 

I do this, the loop is okay i think but why the error is incresing and not decreasing,     

errorrk3:= array(1 .. 7)

 for k from 4 to 10 do

 errorrk3 [k]:=evalf( |(eval(rhs(res), x = 1)) - RK3((x, y) -> -y, 0, 1, 1, k)|)

if errorrRK3[N] <= 0.001 then break; end if;

end do;