@Carl Love You clearly meant to write LinearAlgebra:-Eigenvectors.

You give us an image and that image doesn't even reveal how RK3 is defined.
To get reasonable help you should present the code as text or as an uploaded worksheet.

I do notice that you insist on taylor(abs(Exact-E),h=0); in order to gain insigt into the order of the error.
But recall that  g(h) = O(h^4) (as h ->0 means that there exist delta>0 and C > 0 such that

abs(g(h)) <= C*h^4  for all h satisfying abs(h) < delta.

Thus it is enough to do

taylor( Exact -E,h=0);

@Markiyan Hirnyk Here is a version which doesn't use polar coordinates. Since x' and y' both will be changed when x^2+y^2=1 and since the changes involve x' and yi for both, we need a temporary variable (temp):

res:=dsolve(sys,numeric,events=[[x(t)^2+y(t)^2-1,[temp=diff(x(t),t),diff(x(t),t)=-2*(diff(x(t),t)*x(t)+diff(y(t),t)*y(t))*x(t)+diff(x(t),t),diff(y(t),t)=-2*(temp*x(t)+diff(y(t),t)*y(t))*y(t)+diff(y(t),t)]]],range=0..200);
p1:=plot(1,theta=0..2*Pi,coords=polar,color=blue):
plots:-odeplot(res,[x(t),y(t)],0..200,axes=boxed,axes=none,refine=1):
plots:-display(%,p1);
plots:-odeplot(res,[x(t),y(t)],0..200,axes=boxed,axes=none,frames=500,numpoints=1000):
plots:-display(%,p1);

I copied the line as

Butcher, map(x->if (x=0) then `` else x end if,Butcher);

It executed without an error and returned as it should:

@Markiyan Hirnyk Sorry, I just made an addition above.

@Markiyan Hirnyk Maybe something like this:

PDEtools:-dchange({x(t)=r(t)*cos(theta(t)),y(t)=r(t)*sin(theta(t))},{diff(x(t),t,t)=0,diff(y(t),t,t)=0},[r(t),theta(t)]);
sysP:=solve(%,{diff(r(t),t,t),diff(theta(t),t,t)});
resP:=dsolve(sysP union {r(0)=.2,theta(0)=.5,D(r)(0)=1,D(theta)(0)=4},numeric,events=[[r(t)=1,diff(r(t),t)=-diff(r(t),t)]]);
p1:=plot(1,theta=0..2*Pi,coords=polar,color=blue):
p2:=plots:-odeplot(resP,[r(t)*cos(theta(t)),r(t)*sin(theta(t))],0..40,axes=none,frames=100):
plots:-display(p1,p2);

And with the range argument and maxfun=0:

resP:=dsolve(sysP union {r(0)=.2,theta(0)=.5,D(r)(0)=1,D(theta)(0)=4},numeric,events=[[r(t)=1,diff(r(t),t)=-diff(r(t),t)]],range=0..100,maxfun=0);
plots:-odeplot(resP,[r(t)*cos(theta(t)),r(t)*sin(theta(t))],0..100,refine=1,axes=none);
p1:=plot(1,theta=0..2*Pi,coords=polar,color=blue):
p2:=plots:-odeplot(resP,[r(t)*cos(theta(t)),r(t)*sin(theta(t))],0..100,axes=none,frames=500,numpoints=1000):
plots:-display(p1,p2);

@sarra Sure. But it is evident from taylor(Exact,h=0) that the first two terms are just E. Thus
E-Exact = O(h^2).

@Swhite Using NonlinearFit on your data xlist, ylist, zlist:

f:=x=a*y^b*z^c;

X:=Vector(xlist);
Y:=Vector(ylist);
Z:=Vector(zlist);
A:=Matrix([Y,Z,X],datatype=float);
res:=Statistics:-NonlinearFit(rhs(f),A,[y,z],output=solutionmodule);
res:-Results();

@sarra 

f:=(x,y)->I*y;
res:=dsolve({diff(y(x),x)=f(x,y(x)),y(0)=-I});
MidpointMethod(f, 0, h,-I ,1); #Using the last version I gave
#However it is silly since we just need this part:
#y[1] := y[0]+h*f(x[0],y[0]);
E:=-I+h*f(0,-I);
Exact:=eval(rhs(res),x=h);
taylor(Exact,h=0); #Shows that the error is O(h^2)
plot(abs(E-Exact),h=0..1);
plots:-loglogplot(abs(E-Exact),h=1e-5..1);

@sarra I don't know if this is part of the exercise (the point?) but obviously the value of y[1] is not unimportant.
An obvious change is to define y[1] by an Euler step like this:

MidpointMethod := proc (f, a, b, N) local x, y, n, h;
   y :=Array(0..N);
   x:=Array(0..N);
   h := evalf(b-a)/N;
   x[0] := a;
   y[0] := 1; #Ought to be one of the input values for the procedure
   x[1] := a+h;
   y[1] := y[0]+h*f(x[0],y[0]); #The change
   for n to N-1 do
      x[n+1] := a+(n+1)*h;
      y[n+1] := y[n-1]+2*h*f(x[n], y[n])
   end do;
   [seq([x[n], y[n]], n = 0 .. N)]
end proc;
##########
#Letting y[0] be an input:
MidpointMethod := proc (f, a, b,y0, N) local x, y, n, h;
   y := Array(0 .. N);
   x:=Array(0..N);
   h := evalf(b-a)/N;
   x[0] := a;
   y[0] := y0;
   x[1] := a+h;
   y[1] := y[0]+h*f(x[0],y[0]);
   for n to N-1 do
      x[n+1] := a+(n+1)*h;
      y[n+1] := y[n-1]+2*h*f(x[n], y[n])
   end do;
   [seq([x[n], y[n]], n = 0 .. N)]
end proc;

Then (assuming that i means the imaginary unit??):
res:=dsolve({diff(y(x),x)=I*y(x),y(0)=-I});
#and
MidpointMethod((x,y)->I*y, 0, 1,-I ,12);
#agrees reasonably well.

@sarra You may want to take a look at

http://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods

There seem to be a discrepancy between the initial description and the one in the code:

$$ y[i] = y[i-2] + 2h f(x[i],y[i]);
versus
y[n+1] = y[n-1] +  2h f( x[n], y[n] );

I suppose it is the latter that you mean?

Also slightly confusing is the two different names you use: Midpoint method and Improved Euler method.

@itsme I get this when adding 'allsolutions'
solve(eqs, [omega[n],theta[n]], allsolutions);

Warning, solutions may have been lost
        [[omega[n] = -0.5240921270,

          theta[n] = 1.149798694 + 3.141592654 _Z1],

          [omega[n] = 0., theta[n] = 3.141592654 _Z1], [

          omega[n] = 0.5240921270,

          theta[n] = -1.149798694 + 3.141592654 _Z1], [

          omega[n] = -1.487355972,

          theta[n] = 0.6003934395 + 3.141592654 _Z1]]

# Indeed solutions are lost. There are infinitely many values for omega[n] as is illustrated by the plot I mentioned.
Also notice the comment by Markiyan Hirnyk. If you replace 0.2 by 1/5 you should see the problem. It is not easy.
## Emphasis added later.

@felipe_p To get an overlay of several plots, say p1,p2, p3, you can use

display(p1,p2,p3);
# or
display([p1,p2,p3]);

If instead you use

display(Array([p1,p2,p3]));

you get an array of plots. Your p is an array.

An overlay of plots in your case is not very exciting since each the plot just a part of p[10], but you get the overlay by
plots:-display(seq(p[i],i=0..10));

@Amber Doesn't Carl Loves's version work as you intended? If not, why not?