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

@Markiyan Hirnyk As I recall it plot3d never belonged to the plots package. Thus plots[plot3d] shouldn't work in any Maple version.

@adel-00 So according to your calculation what is the steady state?

@adel-00 Now what do you think is "the steady state" ?

There is a whole curve consisting of "steady states", equilibria, or whatever you like to call them:
rhs~(dsys);
eqs:=evalindets(%,function(identical(t)),f->op(0,f)) =~0;
solve(%,{x,y,z});
#Thus Re(x) = 0, x=I*s with s real.
#Therefore the curve x=I*s,y=-2*s-5,z=8+20/s , with s<>0 and real, consists of equilibria.
#Check: 
eval(eqs,{x=I*s,y=-2*s-5,z=8+20/s});
simplify(%) assuming real;

###None of the equilibrium points can be asymptotically stable. That is just because any neighborhood of one of the equilibria contains another one (indeed infinitely many):
## Stability (but not asymptotical) cannot be ruled out a priori.
##Here is an illustration:
##Start exactly at one of them:
icsE:=eval({x(0)=I*s,y(0)=-2*s-5,z(0)=8+20/s},s=4);
res:=dsolve(dsys union icsE,numeric,maxfun=0);
res(6);
plots:-odeplot(res,[[t,Im(x(t))],[t,y(t)],[t,z(t)]],0..10,legend=[Im(x(t)),y(t),z(t)]);
##Start slightly off that point:
eps:=.001:
icsEeps:=eval({x(0)=I*s+eps,y(0)=-2*s-5,z(0)=8+20/s},s=4);
res2:=dsolve(dsys union icsEeps,numeric,maxfun=0);
res2(6);
# Notice that we need to replace y(t) by Im(y(t)) (or Re(y(t)) ) in the plot:
plots:-odeplot(res2,[[t,Im(x(t))],[t,Im(y(t))],[t,z(t)]],0..50,legend=[Im(x(t)),Im(y(t)),z(t)],size=[1000,400]);

@9009134 I believe that you can trust the numerical result. The method itself is perfectly legitimate (recall the simple (linear) example I gave above).
I would do a simple test for consistency, by which I mean: If I feed RES as an approximate solution to dsolve/numeric/bvp is the solution recovered?
It is:

plots:-odeplot(RES,[[x,1000*f1(x)],[x,1000*f2(x)],[x,f3(x)]],0..1); p1:=%:
#omega = 948.876522424828
sol:=dsolve(dsys3 union {(D@@4)(f3)(0)=-1}, numeric,approxsoln=RES);
plots:-odeplot(sol,[[x,1000*f1(x)],[x,1000*f2(x)],[x,f3(x)]],0..1); p2:=%:
sol(0);
#omega = 948.876522440829
plots:-display(p1,p2); #Indistinguishable

@9009134 Just copy and paste the following code into a fresh worksheet:

restart;
dsys3 := {15*(diff(f3(x), x, x, x, x, x, x))+16*(diff(f3(x), x, x, x, x))+17*(diff(f3(x), x, x, x, x))+18*(diff(f3(x), x, x))+19*(diff(f3(x), x, x))+20*(diff(f3(x), x, x))+21*(diff(f1(x), x, x, x))+22*(diff(f1(x), x))+23*(diff(f1(x), x))+24*(diff(f2(x), x, x))+25*f2(x)+26*f2(x)+27*f3(x)+28*f3(x)+29*f3(x)+30*f3(x)+f3(x)*omega^2 = 0, diff(f1(x), x, x, x, x)+2*(diff(f1(x), x, x))+3*(diff(f2(x), x, x, x))+4*(diff(f2(x), x))+5*(diff(f3(x), x, x, x))+6*(diff(f3(x), x))+7*f1(x)+f1(x)*omega^2 = 0, 8*(diff(f2(x), x, x, x, x))+9*(diff(f2(x), x, x))+10*f2(x)+11*(diff(f1(x), x, x, x))+12*(diff(f1(x), x))+13*(diff(f3(x), x, x))+14*f3(x)+f2(x)*omega^2 = 0, f1(0) = 0, f1(1) = 0, f2(0) = 0, f2(1) = 0, f3(0) = 0, f3(1) = 0, ((D@@1)(f1))(0) = 0, ((D@@1)(f1))(1) = 0, ((D@@1)(f2))(0) = 0, ((D@@1)(f2))(1) = 0, ((D@@1)(f3))(0) = 0, ((D@@1)(f3))(1) = 0, ((D@@2)(f3))(0) = 0, ((D@@2)(f3))(1) = 0};
sys,bcs:=selectremove(has,dsys3,diff);
extra_bcs:={seq(seq(seq((D@@i)(cat(f,j))(a),j=1..2),i=2..3),a=0..1),seq(seq((D@@i)(f3)(a),i=3..5),a=0..1)};
nops(extra_bcs);
for b in extra_bcs do
  try
    print(b=-1);
    res[b]:=dsolve(dsys3 union {b=-1}, numeric,initmesh=1024,approxsoln=[omega=1,f1(x)=x*(1-x),f2(x)=sin(Pi*x),f3(x)=x*(1-x)]);
  catch:
    print(lasterror);
  end try
end do;

eval(res);
RES:=res[(D@@4)(f3)(0)];
RES(0);
plots:-odeplot(RES,[[x,1000*f1(x)],[x,1000*f2(x)],[x,f3(x)]],0..1);
#### From RES(0) you get that omega = 948.876522424828
####Notice that f1 and f2 have been multiplied by 1000 in the plot.

@Preben Alsholm 
Continuing with your new dsys3 this works for one of the extra bcs:

sys,bcs:=selectremove(has,dsys3,diff);
extra_bcs:={seq(seq(seq((D@@i)(cat(f,j))(a),j=1..2),i=2..3),a=0..1),seq(seq((D@@i)(f3)(a),i=3..5),a=0..1)};
nops(extra_bcs);
for b in extra_bcs do
  try
    print(b=-1);
    res[b]:=dsolve(dsys3 union {b=-1}, numeric,initmesh=1024,approxsoln=[omega=1,f1(x)=x*(1-x),f2(x)=sin(Pi*x),f3(x)=x*(1-x)]);
  catch:
    print(lasterror);
  end try
end do;

eval(res);
RES:=res[(D@@4)(f3)(0)];
RES(0);
plots:-odeplot(RES,[[x,1000*f1(x)],[x,1000*f2(x)],[x,f3(x)]],0..1);

@9009134 dsolve/numeric will allow a parameter in the system if provided with an extra condition that will determine that parameter.
I think the idea is that the extra parameter (omega, say) internally is represented by the ode diff(omega(t),t) =0.
The extra condition will of course be put on the unknown functions (in your case f1,f2,f3) not on omega.

I tried the following on your latest dsys3:
extra_bcs:={seq(seq(seq((D@@i)(cat(f,j))(a),j=1..2),i=2..3),a=0..1),seq(seq((D@@i)(f3)(a),i=3..5),a=0..1)};
nops(extra_bcs);
###
for b in extra_bcs do
  try
    print(b=1);
    res[b]:=dsolve(dsys3 union {b=1}, numeric,initmesh=1024);
  catch:
    print(lasterror);
  end try
end do;

As you will see I didn't have success for any of the 14 choices of extra conditions. (Originally I claimed erroneously that there were 8: there are 14 of the kind I talked about). "matrix is singular" was the error message in each case.

@erik10 I upgraded my Windows 8.1 to Windows 10 a couple of days ago. Maple 2015.1 using worksheet mode and 1-D input works without any problems according to my 2 days of experience.

@9009134 This is a totally different system!

You may find inspiration in this simple example, which you can easily solve symbolically, but here it is done numerically:

restart;
ode:=diff(y(x),x,x)+lambda*y(x)=0;
res:=dsolve({ode,y(0)=0,y(1)=0,D(y)(0)=1},numeric);
res(0);
evalf(Pi^2);
plots:-odeplot(res,[x,y(x)],0..1);
res2:=dsolve({ode,y(0)=0,y(1)=0,D(y)(0)=1},numeric,approxsoln=[lambda=40,y(x)=sin(2*Pi*x)]);
res2(0);
plots:-odeplot(res2,[x,y(x)],0..1);

Your problem is to find a nonzero boundary condition that works for your case. You have plenty to try. If you choose the extra condition to be of the form (D@@k)(f_i)(a) = 0 where k is less than the order of f_i and where a is 0 or 1, I think you have 2*(4+4+6) - 14 = 14 different ones. (I corrected an error in this calculation)!!