We must expect that there are other nontrivial solutions corresponding to other values of L.
To find one more I found it easier to use a 2nd order ode satisfied by y(xi).

UU:=(isolate(sys[2],U(xi)));
eval(sys[1],UU);
ode:=collect((lhs-rhs)(%),diff,factor)=0;
odeL:=subs(y(xi)=L*y(xi),convert(ode,D));
resL:=dsolve({odeL,cond,D(y)(1)=1},numeric,approxsoln=[L=1,y(xi)=(1-xi)*(xi-0.8)]);
resL(.9); #Finds L = 1 approximately (as before)
UU; #(the right hand side can be plotted)
plots:-odeplot(resL,[[xi,y(xi)],[xi,rhs(UU)]]);
## Now looking for another value of L near 4:
resL2:=dsolve({odeL,cond,D(y)(1)=1},numeric,approxsoln=[L=4,y(xi)=(1-xi)*(xi-0.8)*(xi-0.9)]);
resL2(1); #L very close to 4
plots:-odeplot(resL2,[[xi,y(xi)],[xi,rhs(UU)]]);

@RomchikM In the code I gave in my answer above you will see the line
plots:-odeplot(res,[[xi,y(xi)],[xi,U(xi)]]);

This shows the graphs of y and U.

########
The problem you are facing is very similar to the following problem.

restart;
ode:=diff(y(x),x,x)+one*y(x)=0;
bcs:= y(0)=0, y(Pi)=0;

dsolve({eval(ode,one=1),bcs}); #infinitely many solutions (as is well known)
dsolve({eval(ode,one=1+1*10^(-8)),bcs}); #No nontrivial solution
res:=dsolve({eval(ode,one=1+1*10^(-8)),bcs},numeric); #Numerical attempt
plots:-odeplot(res,thickness=3);# same result 0
res2:=dsolve({eval(ode,one=1+1*10^(-8)),bcs},numeric,approxsoln=[y(x)=sin(x)]);
plots:-odeplot(res2,thickness=3); #Tiny, but it resembles a*sin(x) for a very small a
resL:=dsolve({ode,bcs,D(y)(0)=1},numeric);
resL(0); #'one' is very close to 1
plots:-odeplot(resL);



So I would suggest using the results I gave in my answer with the adjusted parameter L. Your input equations contain floats and are surely rounded, so why not?

@MTata94 Sure.
But you could also define the procedure q outside p, but then it shouldn't be declared local to p.
If you have a lot of subroutines used in p you may want to define them outside p and then wrap the whole thing in a module where the subroutines could be declared local. That is another story.

Thanks to all.
From the responses I concluded that submitting an SCR is justified. After all, SCR stands for "Software Change Request". It is not only (I hope) a euphemism for bug report.
The request is that dsolve/numeric with compile=true be made to work without a user with no particular technical computer skills having to do anything out of the ordinary.

I searched MaplePrimes for "compile" and found the following discussion

http://www.mapleprimes.com/questions/138276-Win7-And-Compiler

In that a different issue is raised about making Compiler:-Compile work in the Classical interface on 64 bit machines.

One of acer's workarounds redefines 'ssystem' to `dsolve/numeric/ssystem`.
That made me take a look at that procedure on my 64 bit machine with Windows 10 and Maple 2015.1:

I did as follows:

restart;
eval(`dsolve/numeric/ssystem`);
debug(`dsolve/numeric/ssystem`);
dsolve({diff(x(t),t)=x(t),x(0)=1},numeric,compile=true);

I got the following from the last command:

{--> enter dsolve/numeric/ssystem, args = cl -c -DX86_64_WINDOWS -D_M_X64 -we4028 -we4029 -I\"C:\\Program Files\\Maple 2015\\extern\\include\" /Oiyt -DMSVC -Gz -MD -Gy -nologo \"C:\\Users\\Bruger\\AppData\\Local\\Temp\\m2c98262.c\" -F\"oC:\\Users\\Bruger\\AppData\\Local\\Temp\\m2c98262.obj\", 1
                           "ode2.dll"
proc()  ...  end;
                            [1, ""]
<-- exit dsolve/numeric/ssystem (now in comp_proc) = [1, ]}
###
and then the error message reported earlier.

Actually, I don't understand any of it, but maybe somebody else does?
##
I located the file m2c98262.obj in the Temp directory given above. Opened it in Notebook. It is quite readable and starts with:
/**************************************************************
   External C module
   Code translation and wrapper code:
**************************************************************/

@Markiyan Hirnyk Yes, I should have added that on another computer, a 32 bit, on which I also installed Windows 10, I had no problems. This was using Maple 12 though. There Watcom is used.

@Carl Love Sorry, I'm using Maple 2015.1.

Your code for the loop is broken.
I would strongly advise you to use 1-D math input and worksheet interface whenever you want to do anything serious in Maple, i.e. if your sole purpose is not to make things look pretty.

Go to Tools/Options/Display/Input display and pick Maple Notation (that is 1D math notation).
Right after that continue in the Options dialogue box to the Interface tab. Go to Default format for new worksheets.
Choose Worksheet.
Finally press the Apply Globally button at the bottom.
These settings can always be undone if you don't like them.
Notice that changes only show up in new worksheets.

@Boby An iteration method that depends on successively integrating symbolically is in general going to break down very fast because of the inability to complete the integration.
Just take the case g(y) = y*ln(y).
The following will break down already with m=3 even if you replace YF(0) with the limit from the right.

restart;
g:=y->y*ln(y);
Y[0]:=exp(1); #Example chosen for relative simplicity
for m from 1 to 2 do
  YF:=unapply(Y[m-1],x);
  #Y[m]:=limit(YF(x),x=0,right)-int(t^2*(1/t-1/x)*g(YF(t)),t=0..x) assuming x>0,x<sqrt(6)
  Y[m]:=YF(0)-int(t^2*(1/t-1/x)*g(YF(t)),t=0..x) assuming x>0,x<sqrt(6)
end do;
evalf(Y[2]);
collect(%,ln);
#Compare with numerical solution:
ode:=diff(y(x),x,x)+2/x*diff(y(x),x)+y(x)*ln(y(x))=0;
res:=dsolve({ode,y(0)=exp(1),D(y)(0)=0},numeric);
plots:-odeplot(res,[x,y(x)],0..50);
plot(Y[2],x=0..sqrt(6));
plots:-odeplot(res,[x,y(x)-Y[2]],0..sqrt(6),thickness=3);

@adel-00 You obviously are not considering the real version of the original complex. The first equation in the original complex version had epsilon added, here it is multiplied (for one difference).

@adel-00 As I said earlier (for the case epsilon=5,Delta1=2,Delta2=4) the stady states are
x=I*s,y=-2*s-5,z=8+20/s, where s<>0.
Your version written as

X∞=− iε∆1/∆1∆2+c

Y∞ = -εc / ∆1∆1 +c

Z∞ = c 

I don't quite understand. Could insert multiplication signs and parentheses as needed?
##Added: Changing my own parametrization (in s) to parametrization in z I get
[x = -I*Delta2*epsilon/(Delta1*Delta2+z), y = -epsilon*z/(Delta1*Delta2+z), z = z]
####

As I said none of those can be asymptotically stable. Thus you will never get close to reaching any one of them by integrating for even a long time (unless you start at the point itself). (see modification below).
A closer examination will show (I believe) that none of the points are even stable.

## Note added:
Instead of saying that solutions should approach "the steady state" you may have meant that the solutions should approach the set of steady states?
As in the following very simple example:
sys:={diff(x(t),t) = -x(t), diff(y(t),t) = 0};
The set of equilibria is the whole y-axis. Each solution  (x(t),y(t) tends to (0,y0) for some y0 (y0 dependent on the given solution!)
Whether that kind of behavior occurs in your system is at least questionable.
Do you have any reason to believe that any solution will approach some steady state (which may vary according to the given solution)?

The problem of isolating the highest derivatives can also be solved directly by differentiating one of the equations yourself as in this approach.

restart;
odes:=diff(x1(t),t)*diff(x2(t),t$2)*sin(x1(t)*x2(t))+5*diff(x1(t),t$2)*diff(x2(t),t)*cos(x1(t)^2)+t^2*x1(t)*x2(t)^2=exp(-x2(t)^2),diff(x1(t),t$2)*x2(t)+diff(x2(t),t$2)*diff(x1(t),t)*sin(x1(t)^2)+cos(diff(x2(t),t$2)*x2(t))=sin(t);
ics:=x1(0)=1,D(x1)(0)=1,x2(0)=2,D(x2)(0)=2;

ode1:=isolate(odes[2],diff(x1(t),t,t));
ode2temp:=subs(ode1,odes[1]);
diff(ode2temp,t);
ode2:=isolate(%,diff(x2(t),t,t,t));
#We need an initial condition for(D@@2)(x2):
convert(ode2temp,D);
ic2:=(lhs-rhs)(eval[recurse](%,{t=0,ics}));
plot(subs((D@@2)(x2)(0)=yp4,ic2),yp4=-5..5); #3 solutions
ics2List:=(D@@2)(x2)(0)=~sort([solve(evalf(ic2),(D@@2)(x2)(0))]); #3 possible initial conditions
res:=seq(dsolve({ode1,ode2,ics,ics2List[i]},numeric),i=1..3); #All 3 at once as a sequence
plots:-odeplot(res[1],[seq([t,diff(x1(t),[t$i])],i=0..1)],0..0.212);
plots:-odeplot(res[1],[seq([t,diff(x2(t),[t$i])],i=0..2)],0..0.212);
plots:-odeplot(res[2],[seq([t,diff(x1(t),[t$i])],i=0..1)],0..0.218);
plots:-odeplot(res[2],[seq([t,diff(x2(t),[t$i])],i=0..2)],0..0.218);
plots:-odeplot(res[3],[seq([t,diff(x1(t),[t$i])],i=0..1)],0..0.252);
plots:-odeplot(res[3],[seq([t,diff(x2(t),[t$i])],i=0..2)],0..0.251);

@Al86 To plot the actual error you would have to know the exact solution, but you don't. To get an idea of the success of this you could try doubling n as is done below and then plot the result Ymono2-Ymono:

n:=20:
Y:=t->1+add(c[i]*t^i,i=1..n);
eqt:=eval(eq2,y=Y):
T:=[seq(0.05..1,0.05)];
nops(T); #Same as n
eqs:={seq(evalf(eval(IntegrationTools:-Expand(eqt),t=t0)),t0=T)}:
sol:=fsolve(eqs);
eval((rhs-lhs)~(eqs),sol); #Checking collocation
Ymono2:=eval(Y(t),sol); #The solution using the monomials t^i
plot(Ymono2,t=0..1,size=[1800,400]);

plot(Ymono2-Ymono,t=0..1,size=[1800,400]);

@Boby I responded to that below regarding your Lane-Emden equation.

@Boby In my response to your question in
http://www.mapleprimes.com/questions/204964-ODE-Of-LaneEmden-Type

I gave a simple version for the equation diff(y(x),x,x)+2/x*diff(y(x),x)+exp(-y(x)) = 0 with y(0)=0.

restart;
Y[0]:=0: #Numbers can work as constant functions
n:=3:
for m from 1 to 3 do
  YF:=unapply(Y[m-1],x);
  Y[m]:=-int(t^2*(1/t-1/x)*exp(-YF(t)),t=0..x)
end do;

The problem we run into there is that pretty soon (for Y[3]) the integrals cannot be done symbolically. That is no problem for your present ode, which is linear, homogeneous, and has constant coefficients. It doesn't get much simpler.