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;

@k_1123 Certainly. An easy implementation (but not necessarily the most efficient) would be:

f:='f':
Student:-Calculus1:-ApproximateInt(f(x), a..b, method = simpson[3/8],partition=1);
F:=unapply(subs(b=a+3*h,%),a,h);
# "Stealing" the n=1 formula from Maple.

simp:=proc(a,b,n) local h,i,F,x0,S;
  F := (a,h)->(3/8)*h*(f(a)+3*f(a+h)+3*f(a+2*h)+f(a+3*h)); #n=1 formula
  h:=(b-a)/n/3;
  x0:=a;
  S:=F(x0,h);
  for i from 2 to n do
    x0:=x0+3*h;
    S:=S+F(x0,h) 
  end do;
  S
end proc;

simp(a,b,1);
simp(a,b,2);
factor(%);

I think you are using the wrong formula. You can see a single step (or more) in the abstract by doing:

f:='f':
Student:-Calculus1:-ApproximateInt(f(x), a..b, method = simpson[3/8],partition=1);
Student:-Calculus1:-ApproximateInt(f(x), a..b, method = simpson[3/8],partition=2);
Student:-Calculus1:-ApproximateInt(f(x), a..b, method = simpson[3/8],partition=3);
Clearly your code does not do this job.

@teolog Given the opportunity I cleaned up the code some and I'm back to using U (instead of P).
I give 2 versions of using events.
I also show how to use the parameters option in dsolve. That is particularly useful when you want to examine dependence on parameters, in your case the initial values. I show an example of an animation in U0 holding x10 and x20 fixed.

restart;
params:=[x10=0.2, x20=0.4, U0=0.25];
eq3a :=(1/4)*(U(t)^2*x1(t)-2*U(t)*x1(t)*sqrt(x2(t))-U(t)^2*x2(t)+2*U(t)*x2(t)*sqrt(x1(t))-2*(diff(U(t), t))*x1(t)^(3/2)*U(t)+4*(diff(U(t), t))*x1(t)^(3/2)*sqrt(x2(t))+2*(diff(U(t), t))*x2(t)^(3/2)*U(t)-4*(diff(U(t), t))*x2(t)^(3/2)*sqrt(x1(t)))/(x1(t)^(3/2)*x2(t)^(3/2));
sysa := diff(x1(t), t) = U(t)-sqrt(x1(t)), diff(x2(t), t) = U(t)-sqrt(x2(t));
eq3:=isolate(eq3a=0,diff(U(t),t));
sys:= eq3,sysa;
res:=dsolve(eval({sys,x1(0)=x10,x2(0)=x20,U(0)=U0},params),numeric);
plots:-odeplot(res,[[t,U(t)],[t,x1(t)],[t,x2(t)]],-3..8,view=0..1);
resE1:=dsolve(eval({sys,x1(0)=x10,x2(0)=x20,U(0)=U0},params),numeric,events=[[U(t)-0.35,halt]]);
resE1(4);
plots:-odeplot(resE1,[[t,U(t)],[t,x1(t)],[t,x2(t)]],-3..4,view=0..1);
resE2:=dsolve(eval({sys,x1(0)=x10,x2(0)=x20,U(0)=U0,b(0)=0},params),numeric,discrete_variables=[b(t)::float],events=[[U(t)-0.35,b(t)=t]]);
resE2(4);
plots:-odeplot(resE2,[[t,U(t)],[t,x1(t)],[t,x2(t)],[t,b(t)]],-3..4,view=0..1.3,thickness=3,legend=[U,x1,x2,b]);
resE2(2);
#####################################################
respar:=dsolve({sys,x1(0)=x10,x2(0)=x20,U(0)=U0},numeric,parameters=[x10,x20,U0]);
respar(parameters=params);
plots:-odeplot(respar,[[t,U(t)],[t,x1(t)],[t,x2(t)]],-3..4,view=0..1);
Q:=proc(x10,x20,U0) respar(parameters=[x10,x20,U0]); respar end proc;
plots:-animate(plots:-odeplot,['Q'(0.2,0.4,U0),[[t,U(t)],[t,x1(t)],[t,x2(t)]],0..8],U0=0.1...0.5);

@teolog You can use events to stop integration at say U(t) =0.35: Just do (remember that I called U for P):

res:=dsolve(NewSys union {x1(0)=x10,x2(0)=x20,P(0)=U0},numeric,events=[[P(t)-0.35,halt]]);

It does seem that the solution P(t() appoaches a limit though, as t -> infinity.

@Markiyan Hirnyk U (or as I call it P) is determined by U(0) = U0 in any case.

I really don't understand what you are asking. Can you put it in a different way?

@Carl Love Just a comment:
Surely, it is documented that discont should be either true or false (or a list).
However, plot would not have accepted discont=FAIL if Plot:-Preprocess had used type {truefalse,list} instead of {boolean,list} as the type check for discont.

plot(1/x, x=-1..1, discont= xxx ,thickness=2, color = blue, legend = '1/x');
Error, (in plot) invalid input: Plot:-Preprocess expects value for keyword parameter discont to be of type {boolean, list}, but received xxx
#So the following goes through:
Plot:-Preprocess(1/x, x=-1..1, discont= FAIL ,thickness=2, color = blue, legend = '1/x');
#The procedure Plot:-Preprocess:
showstat(Plot:-Preprocess);

In Maple 17.02 I get:

restart;
int(1/(sin(t)+1/2), t = -1/2*Pi .. 1/2*Pi);
                           undefined

i.e. as it should be. The same in Maple 16.02.
In Maple 15 though, the answer is wrong.

@Andriy What I was saying was that when using the selection operation [] (help page: ?selection) in our case on A, A is not fully evaluated first and then the selection made.

When i = 2 the string of events appears to be A[2]  ->  {a,b,c}[2] -> b -> a.
Normally Maple evaluates fully, but there are important exceptions, think of procedures or tables.
Try evaluating A just once, i.e. not fully using eval(A,1):

restart;
A := {a, b, c};
b := a;
A;

eval(A,1);
eval(A,2);
#Here the normal full evaluation is used:
seq(op(i,A),i=1..3); #Error because of nops(A)=2
seq(op(i,A),i=1..2);

########### Another example:
restart;
A := [a, b, c];  #or try A:=a,b,c; instead
b := NULL;
A;

eval(A,1);
eval(A,2);
seq(A[i],i=1..3);