@siamak taghavi Numerical solution clearly requires sigma01[xx] and sigma01[xy] to be numerical. In your worksheet you substitute values for these in your expression for the symbolic solution. Those are the ones I used in the numerical solutions res1 and res2. So I fail to understand why you say that is wrong.

@siamak taghavi Your worksheet contains no text, i.e. no explanation whatsoever. Thus it takes time to find out what is going on and what you think is the problem.
As far as I understand it it, you are solving two boundary value problems on [-a,a] for the two odes. 
b11 and b21 are supposed to be phi'''(rho) and psi'(rho) for the first set of boundary values, whereas b12 and b22 are the corresponding values for the second set.
Your immediate problem in using a numerical method is that rho > a, thus outside [-a,a].
Therefore you must first solve the bvp-problem, find the relevant function values at the point xi=a, and then solve an initial value problem starting at xi=a.
Doing that, I get:

res1:=dsolve({VR2=0,VS2=0,phi(a)=1,phi(-a)=1,D(phi)(a)=0,D(phi)(-a)=0,psi(a)=0,psi(-a)=0},numeric);
plots:-odeplot(res1,[xi,phi(xi)],-a..rho); #Notice rho and the complaint.
res1(rho);
res1(a);
phi2a,phi3a,psi1a:=op(subs(res1(a),[diff(phi(xi),xi,xi),diff(phi(xi),xi,xi,xi),diff(psi(xi),xi)]));
ics1a:={phi(a)=1,D(phi)(a)=0,(D@@2)(phi)(a)=phi2a,(D@@3)(phi)(a)=phi3a,psi(a)=0,D(psi)(a)=psi1a};
res1cont:=dsolve({VR2=0,VS2=0} union ics1a,numeric);
res1cont(rho);
b11,b21:=op(subs(res1cont(rho),[diff(phi(xi),xi,xi,xi),diff(psi(xi),xi)]));

@OtissKiller For convenience I repeat the code, but I have made some additions:

restart;
DE1 := diff(Y(t), t) = 5*Y(t)*ln(b(t)/Y(t))-5*Y(t);
DE2 := diff(b(t), t) = 5*b(t)*Y(t)^(3/2)-5*Y(t);
res:=dsolve({DE1,DE2,Y(0)=1,b(0)=2},numeric);
plots:-odeplot(res,[[t,b(t)],[t,Y(t)]],-.3..0.5);
plots:-odeplot(res,[Y(t),b(t)],-.3..0.5);
solve(rhs~({DE1,DE2}),{b(t),Y(t)});
allvalues(%);
pt:=subs(%,[Y(t),b(t)]); #The equilibrium point
evalf(%);
#A phaseplot with orbits starting now at 10 different points.
DEtools[DEplot]({DE1,DE2},[Y(t),b(t)],t=-1..10,
[
[Y(0)=0.5,b(0)=1.4],
[Y(0)=0.52,b(0)=1.4],
seq([Y(0)=1,b(0)=.5*i],i=1..4),
seq([Y(0)=.1,b(0)=.4*i],i=1..4)
],
Y=0..2,b=0..2,stepsize=0.01,linecolor=blue); p1:=%:
#A plot of the equilibrium point:
p2:=plot([pt],style=point,symbolsize=30,symbol=solidcircle,color=black):
plots:-display(p1,p2); #Displaying the two plots together
#What type of point is pt? From the phase plot it looks like a saddle point:
RHS:=subs(Y(t)=Y,b(t)=b,rhs~([DE1,DE2]));
J1:=VectorCalculus:-Jacobian(RHS,[Y,b]); #The jacobian
J:=unapply(J1,Y,b): #The jacobian as a function of Y and b
J(op(pt)); #The jacobian at pt
LinearAlgebra:-Determinant(%);
#The determinant is negativ, thus one eigenvalue is negative, the other positive. The point pt is indeed a saddle point.

@radzys You should use PHi1, which is the numeric procedure for computing the first derivative of phi.
PHI is the numerical procedure for computing phi. If you want to use PHI for computing the derivative of phi at 20, you should use fdiff. But it is a bad idea, since you already have it in PHI1(20).
#Results:
PHI1(20);
                  -1.5952038036214675
fdiff(PHI(x),x=20);
                          -1.595203815

@Sukhdarshan You can modify PC to return unevaluated if the input is not numeric. That is helpful on occasion.
PC:=proc(s,aa) local k,r,N;
   if not type([s,aa],list(numeric)) then return 'procname(_passed)' end if;
   N:=5;
   r[0]:=0;
   for k from 1 to N do
      r[k]:=RootFinding:-NextZero(unapply(eval(h,{S=s,a=aa}),P),r[k-1]);
      if r[k]=FAIL then N:=k-1; break end if
   end do;
   [seq(r[i],i=1..N)]
end proc;
#Now you can try
plot(PC(0.5,a)[1],a=0.01..0.99);q1:=%:
plot(PC(0.5,a)[2],a=0.01..0.99);q2:=%:
plots:-display(q1,q2);
#As for the rest of your worksheet you can do something like this:
W11,W22,W13,W23:=op(map(eval,[V11,V22,V13,V23],{a=aa,S=s,P=PC(s,aa)[2]}));
Q1:=-Q3*W13/W11:
Q2:=-Q3*W23/W22:
W1,W2:=op(eval([w1,w2],S=s)):
#For s=0.5 and aa=0.2 you get:
plot(eval(W1,{Q3=1,s=0.5,aa=0.2}), x1 = 0 .. 0.5, color = red); p1:=%:
plot(eval(W2,{Q3=1,s=0.5,aa=0.2}), x2 = 0 .. 0.5, color = blue); p2:=%:
#which you can display together as you did.
############
As for getting an anytical answer for the solution of h=0 for P in terms of a (and S) you won't get very far. The best you can do is to get a RootOf answer, which is not useless, but for plotting purposes certainly implicitplot is way better:
res:=solve(h=0,P);
plot(eval(res,S=0.5),a=0..1); #Fine, but only one branch
#Now inserting the P-interval 0.6..3 into RootOf as a second argument: 
res2:=evalindets(res,specfunc(anything,RootOf),b->subsop(1=(op(1,b),.6..4),b));
plot(eval(res2,S=0.5),a=0..1,adaptive=false); #Not wonderful!

@radzys No, you are not going to get an analytical expression for phi(x), but you have a numerical procedure which will give you phi and its derivatives. Just use output=listprocedure:
res1a:=dsolve(eval(dsys,{k_r=k_h/2}) union {phi(0)=100},numeric,output=listprocedure);
PHI,PHI1,PHI2,PHI3:=op(subs(res1a,[seq(diff(phi(x),[x$i]),i=0..3)]));
plot(PHI,0..88);
PHI(81.12345);
PHI(x); #returns unevaluted as it should.
##Using the previous res1 and attempting to solve symbolically (even without boundary conditions):
ode:=eval(EOM,res1(0)[-1]);
dsolve(ode);
#The result is a DESol structure, see the help page
?DESol

@siamak taghavi You could just upload relevant files here on MaplePrimes.

@radzys You didn't give us the value of k_r.
Presumably there are several (infinitely many) eigenvalues?

@Carl Love Since the system for given values of omega and k_r clearly has the trivial solution phi(x)=0 (all x), we may assume that the problem is intended as an eigenvalue problem, i.e. omega has to be determined too. We still need to know the value of k_r.
With the syntax corrected (mainly EI and A defined as functions), and with omega^2 replaced by lambda, I tried with the extra condition phi(0)=100 (to determine lambda):

res:=dsolve(eval(dsys,{k_r=k_h/2}) union {phi(0)=100},numeric);
which gave me lambda = 1.10074762511253*10^(-10).

@siamak taghavi If I understand your data in the file file-1.xlsx correctly then the first column consists of x-values and the two other columns are the corresponding phi and psi values.
Why do the x-values continue past the interval -1.366..1.366?
(and also why are there no negative values of x?)
You are talking about a second solution. Does this mean that this problem has multiple solutions? Since the system is linear this would imply infinitely many!

Calling the sequence of ODEs sys I tried setting all parameters to 1:
params:=indets({sys},name) minus {t} =~1;
eval({sys},params);
dsolve(%);
In view of the complexity of the output I think it doubtful that you could get a symbolic answer for your original problem.

This is really a continuation of
http://www.mapleprimes.com/questions/201770-How-To-Override-Plus--Operator-So-That

@rit I'm not sure I understand what you want to accomplish. Certainly in using matrix multiplication addition is performed, but that would (luckily) not be affected by the device I gave. If it did quite a lot of stuff in Maple would similarly be affected.
A simple example:
restart;
local `+`;
`+`:=proc(a,b) :-`+`(a^~2,b^~2) end proc;
mtaylor(exp(x),x=0,2);
                             1 + x
%;
                             1 + x
eval(%,x=2);
                               3
1+x;
                              2    
                             x  + 1
eval(%,x=2);
                               5

Basically invfourier cannot do the job. invfourier calls fourier and that cannot do it. That it returns 0 is not optimal of course.
I tried replacing all floats by 1:
u:=evalindets(vout_fourier_num2,float,1);
#Then you are dealing with
#u := fourier(phi[3](t), t, w) = w^2*cosh(w)*exp(1+w*(I+w))/(w^4+I*w^3+w^2+I*w+1)
#If you try
invfourier(u,w,t);
# or
fourier(u,w,t);
# you will see that they return partially unevaluated.

@Carl Love I tried some (unsophisticated) code I had written recently. It returned a solution which didn't resemble the one obtained by continuation in lambda. That made me suspicious. Then I recalled reading about the equation of Blasius (and Falkner–Skan) not having unique solutions for certain values of the parameter (Philip Hartman: ODEs, Harold Davis, Intro. to Nonlinear Diff. and Int. Equations, Wikipedia). The present equation is from the same area. 
Hartman imposes the extra condition that the derivative be positive and bounded (by 1) for all t (~eta). A similar condition here would weed out some solutions, but even then it is not clear to me if there is only one. Obviously in Hartman inf = infinity, not 10.

An extra problem here is that there is an apparent singularity at values for theta where 1 - f(eta)^2*lambda = 0, since that quantity appears as the denominator when isolating for f'''.