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

@J4James With z=P/L^2*exp(-L*x), T(x) = U(z)  and bcs:=T(0)=1,T(infinity)=0;

we have that x = 0 corresponds to z = P/L^2.
But x = infinity corresponds to z = 0 only if L > 0.
If L < 0 then x = infinity corresponds to z = infinity, since then exp(-L*x) -> infinity as x -> infinity.
(I'm assuming P > 0 in this the latter case too. Otherwise x = infinity corresponds to -infinity).

@J4James I think you miscalculated some.
But the following does come up with a relatively compact answer. The boundary conditions for U are assuming that L > 0.

restart;
ODE:=(diff(T(x), x, x))+P*(S+a*(1-exp(-L*x))/L)*(diff(T(x), x))=0;
bcs:=T(0)=1,T(infinit)=0;
var:=z=P/L^2*exp(-L*x);
var2:=x=solve(var,x);
PDEtools:-dchange({var2,T(x)=U(z)},ODE,[z,U(z)]);
collect(%,diff,factor);
ode1:=isolate(%,diff(U(z),z,z));
dsolve({ode1,U(P/L^2)=1,U(0)=0});

It seems to me from the image you presented that the first x is right next to the first parenthesis (.
This means that there is a missing multiplication sign. The same appears to be true for the first y.

@J4James Yes indeed. The exponential exp(4.77*z) seems to be the underlying problem.

Try with infinity = 1 just to see what is going on.
I noticed that your L is in fact negative, so I also tried replacing L with -L1 in in analytical computation of the solution. The solutions returned in the two cases used different special functions. THe one with Whittaker performed better at the end.

restart:
ODE:=diff(T(z),z$2)+A1*(S-1/L+1/L*exp(-L*z))*diff(T(z),z)+A2*T(z)=0;
bcs:=T(0)=1,T(inf)=0;
ODE2:=subs(L=-L1,ODE);
sol2:=dsolve({ODE2,bcs});
sol1:=dsolve({ODE,bcs});
################
L:=(S-sqrt(S^2+4*alpha*epsilon))/(2*epsilon);
A1:=3*R*P/(epsilon2*(3*R+4));
A2:=3*n*R*P/(epsilon2*(3*R+4));
L1:=-L;
plot(eval(rhs(sol2),{epsilon=.01,epsilon2=.1,alpha=5,R=1,n=1,P=6,S=1,inf=1}),z=0..1);
plot(eval(rhs(sol1),{epsilon=.01,epsilon2=.1,alpha=5,R=1,n=1,P=6,S=1,inf=1}),z=0..1);
########Numerically
Digits:=15:
ODE;
BVP:=eval({ODE,bcs},{epsilon=.01,epsilon2=.1,alpha=5,R=1,n=1,P=6,S=1,inf=1});
sol:=dsolve(BVP,numeric,maxmesh=1024);
plots:-odeplot(sol,[z,T(z)],0..1);

@Luka You got infinitely many solutions with 4 free parameters, although one of them t4 must be positive as we noticed. With res being defined as above you can find the free parameters as the ones appearing in res as ti = ti:

lhs~(select(evalb,res));

I get {t3, t4, t6, t9}.

@J4James Well, the main problem is the unevaluated limit in the result. That is clearly due to the boundary condition at infinity. Replace it (as in the original question) by a symbol, say inf, then

bcs:=T(0)=1,T(inf)=0;

sol:=dsolve({ODE,bcs}):
plot(eval(rhs(sol),{epsilon=.01,epsilon2=.1,inf=5}),eta=0..5);

Whether or not the result contains imaginary numbers is irrelevant. Think of
I*(exp(I*t)-exp(-I*t)) for t real.