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

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