@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?

@teolog Using the notation from the full code I posted you could do:

eq3p:=diff(U(t),t)=piecewise(U(t)<Umax,rhs(eq3),0);
params2:={op(params),Umax=0.35};
sysP:=eq3p,sysa;
resP:=dsolve(eval({sysP,x1(0)=x10,x2(0)=x20,U(0)=U0},params2),numeric);
plots:-odeplot(resP,[[t,U(t)],[t,x1(t)],[t,x2(t)]],-3..4,view=0..1);

@henrylyl You are assigning to the formal parameters X and Y.

Try (with your global X):
p:=proc(Z) Z:=<0,0,0,1,1,1> end proc;
p(X);
p(aaa); #aaa just a name
aaa;
p2:=proc(n) n:=5 end proc;
p2(7);
p2(w); #w a name
w;



The problem is easily correceted:

Make these changes (only):

SLRrepeatedsample:=proc(X1,Y1,N,CI)
local x, y, n, ymu, ySE, yvar, x2,b,bCI,i,X,Y;
X:=X1;
Y:=Y1;