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

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