Your system is linear in x and y, so certainly you can trust solve to have found all solutions of this simple system. What does t have to do with this system?

@abscissa If you have several equilibrium points for a given set of concrete values of the parameters (i.e. in your case A, B, C, D, E, F, G, H, S, T, U, V, W, Y, x) then you may find that defining J as a function of the variables a, b, c, d, e, f, g, h, s, t, u, v may be convenient.
In two steps for clarity:
J1 := VectorCalculus:-Jacobian([eq_1, eq_2, eq_3, eq_4, eq_5, eq_6, eq_7, eq_8, eq_9, eq_10, eq_11, eq_12], [a, b, c, d, e, f, g, h, s, t, u, v]);
J:=unapply(J1,[a, b, c, d, e, f, g, h, s, t, u, v]):
#Now you can use J as a function as in
pt:=seq(1..12);
J(pt);

@pedromneto If there is a way to solve this system by using the laplace transform, I would like to hear about it!

@RemonA What is the capital X appearing in the error message as an argument to sin and cos?

We need more detail than that to give a useful answer.

@roussea18u 
I notice several problems among which are:
1. F is differentiated w.r.t. the variable in the form diff(F,t). This makes sense only if F is an expression in t, not a function. But you assign to F(0) which means that F is supposed to be a function. Also what is F?
Presumably it is globally given and is the right hand side of the ode? Or is it the F appearing on page 10 of the paper in the link. In that case it is a function of two variables?
2. Now if t is time and is used as above, then you cannot assign t to be a vector. Also that assignment is strange. What is t:=vector(i,1) supposed to do? i is unassigned at that point. (Actually you have =, not :=, but then nothing is being assigned to, but I assume that was the intention. There are more cases  of equality signs which should be assignments).
3. You need Newmark:=proc(...) local ... : Notice := and no semicolon.

My suggestion actually is to give us the differential equation including initial conditions. The we can let Maple's dsolve do the problem.

You provide an image. It would be so much easier for us to analyze your code if you gave us the code as text here in MaplePrimes or if you uploaded a worksheet with the code.

A somewhat less trivial example (Volterra-Lotka model):
sys:=diff(x(t),t)=x(t)*(1-y(t)),diff(y(t),t)=-y(t)*(1-x(t));
res:=dsolve({sys,x(0)=1,y(0)=1/2},numeric,output=Array([seq(.1*i,i=0..100)]));
A:=res[2,1];
plot(A[..,2..3],caption="Plot in phase space",labels=[x,y]);
plot([A[..,[1,2]],A[..,[1,3]]],labels=[t,"x,y"]);
ExportMatrix("F:/testMatrix.txt",A);

@Mazaya Jamil For each concrete value of phi you can use fsolve:

f9:=RootOf(126*exp(I*phi)*_Z^5+(15-351*exp(I*phi))*_Z^4+(107+740*exp(I*phi))*_Z^3+(354-1110*exp(I*phi))*_Z^2+(1056*exp(I*phi)+624)*_Z+480-480*exp(I*phi));
p:=op(f9); #The polynomial in _Z
fsolve(eval(p,phi=.12345),_Z); #Example
#Here is a procedure which computes the roots for each given value of phi:
Q:=proc(phi1) if not phi1::numeric then return 'procname(_passed)' end if;
   fsolve(eval(p,phi=phi1),_Z)
end proc;
Q(8);
Q(.12345);
Q(x);
#It returns unevaluated if it receives non numeric input.
#Notice that in your case you don't need the extra argument 'complex' because the coefficients of p are imaginary.
Here the extra argument is necessary (if you want imaginary roots too):
fsolve(z^2+z+1,z,complex);
############
Plot of all 5 roots as phi varies from 0 to 2*Pi:
For efficiency add 'option remember' to Q:
Q:=proc(phi1) option remember; if not phi1::numeric then return 'procname(_passed)' end if;
   fsolve(eval(p,phi=phi1),_Z)
end proc;
plots:-complexplot([seq(Q(phi)[i],i=1..5)],phi=0..2*Pi,style=point);
#And an animation of the same plot:
plots:-animate(plots:-complexplot,[[seq(Q(phi)[i],i=1..5)],phi=0..pp,style=point,numpoints=round(pp*5+1)],pp=0..2*Pi);

Try leaving out the initial condition:
pdsolve(sys[1]);

If that is correct it is no wonder that the initial value problem cannot be solved.

Not knowing very much about MatLab I was wondering if Maple's 'break' is what you need:

?break

In the help page you find the example

L := [1, 2, "abc", "a", 7.0, infinity]:
for x in L do
    if type(x, 'string') then
        print(x);
        break;
    end if;
end do;

@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