In order to help we need details: the odes and the initial conditions as text or as an uploaded worksheet.

@J4James The example you give is an example for which no point computed is real. Therefore that warning and an empty plot. If on some interval lambda1 is real then plot does give a graph on that interval and nothing elsewhere: see my example above.
But if for some reason you prefer not to get the warning you could do
cdn:=alpha>=-(1/2)*S^2/(K+2);
eval(cdn,{alpha=-1,S=-0.5}); #Never satisfied for K=0..2
plot(eval(piecewise(cdn,lambda1,undefined),{alpha=-1,S=-0.5}),K=0.0..2);
but I fail to see the reason why.

@Markiyan Hirnyk I have seen you use that phrase "empty words" quite a few times. I'm not sure whether I should be offended or not. Maybe some time at your leisure you could explain what you mean by that.

But to prove the converse:
if lambda1 is a root of order >= k+1 in the polynomial p(lambda) then
p(lambda) = q(lambda)*(lambda-lambda1)^(k+1) where q(lamda) is a polynomial.
Thus (D@@i)(p)(lambda1) = 0 for i=0..k. Thus from "that equation" i.e.
 "p(D)(t^k*exp(lambda*t) ) = kth derivative of p(lambda)*(exp(lambda*t) wrt. lambda."
follows that

p(D)(t^k*exp(lambda1*t) ) = 0 for all t.

You say "My last question is: "  and yet you don't formulate a question, but repeat the assignment you have been given. I assume that you are not expecting us to do your assignment.

@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