thanks for your help and support,but there is a little problem,and that is that when i use this procedure to solve my own problem,or at this problem,the matrix K1 is also fnction of b[i,j]s, but i want them to be solved just as functions of a[i]s,becuase i will assign nemerics to this a[i] and want to use them in a loop,so that i will need the inverse of the matrix K1,which in a 18*18 matrix take so many time if all of the elements are not numeric
what should i do? could u plz help?thanks alot,and here is my own code : solve,K.mws

@Doug Meade no,i think  we can write them as MATRIX multiply to a VECTOR but i do not know how?, the matrix itself can be the function of the coefficients,thnx anyway.

@Markiyan Hirnykbut i thought u have took the file,and i had to do that,thnx anyway,

@Markiyan Hirnykbut i thought u have took the file,and i had to do that,thnx anyway,

@Markiyan Hirnyk i can not solve my equations since when using dsolve,it says that it has singularity,and when i want to solve them,using stiff=true options and maxfun=0,it takes more than a day to solve it.what can i do to solve them correctly?do u see any problem in the form of my equations?

when i see i can not solve it using dsolve,i want to solve it using step by step integration,and i need to convert these equations to vector mode,but now if u see there is a problem,all of components of the first equation are not included in the matrix u created,for example in the element first equations,all of multipleres of a[13] are not included,if u see in the pic,the coeeficient 0.1535e12*a[13]^2 is not in the line 13 of the matrix A(13,1)? i hope it will be clear enough,thnx.

@Markiyan Hirnyk i can not solve my equations since when using dsolve,it says that it has singularity,and when i want to solve them,using stiff=true options and maxfun=0,it takes more than a day to solve it.what can i do to solve them correctly?do u see any problem in the form of my equations?

when i see i can not solve it using dsolve,i want to solve it using step by step integration,and i need to convert these equations to vector mode,but now if u see there is a problem,all of components of the first equation are not included in the matrix u created,for example in the element first equations,all of multipleres of a[13] are not included,if u see in the pic,the coeeficient 0.1535e12*a[13]^2 is not in the line 13 of the matrix A(13,1)? i hope it will be clear enough,thnx.

@Markiyan Hirnyk could u plz email me so that i can send the file to your email?thank u for your attention ,

@Markiyan Hirnyk could u plz email me so that i can send the file to your email?thank u for your attention ,

could u  rewrite the code after edition?i could not understand what u mean by this edit?,
this code has another problem also, in matrix A  some expressions are not included,the expression including the power of 2 of the function,for example :
eq1:=.1689228391e12*a[13][t]^2+6406045.412*a[14][t]^2-4317791.317*a[15][t]^2+2291210.983*a[16][t]^2-.2842690977e11*a[17][t]^2-.1782456232e12*a[18][t]^2+.9846526429e12*a[1][t]+.2533881291e12*a[2][t]+.3076607771e11*a[3][t]+8105875.203*a[4][t]-5054889.363*a[5][t]-34561.30764*a[6][t]+9251036.956*a[15][t]*a[13][t]-19358824.20*a[14][t]*a[15][t]+.4614468816e11*a[13][t]*a[14][t]-.1082338049e11*a[18][t]*a[15][t]+23266721.99*a[18][t]*a[13][t]+25978444.80*a[18][t]*a[14][t]-.5529464441e11*a[18][t]*a[16][t]-.1615478137e12*a[17][t]*a[18][t]-27435155.86*a[16][t]*a[15][t]-6275707.162*a[16][t]*a[13][t]-20873274.82*a[14][t]*a[16][t]+37522329.49*a[15][t]*a[17][t]-5539558.102*a[17][t]*a[13][t]+26923086.46*a[14][t]*a[17][t]-.1158227844e11*a[16][t]*a[17][t]:

 but the terms of A(1,i) i=1..18 are :

.9846526429e12
.2533881291e12
.3076607771e11
0.8105875203e7
-.5054889363e7
-34561.30764
0
0
0
0
0
0
23266721.99*a[18][t]-6275707.162*a[16][t]+.4614468816e11*a[14][t]-5539558.102*a[17][t]+9251036.956*a[15][t]
.4614468816e11*a[13][t]-20873274.82*a[16][t]+25978444.80*a[18][t]-19358824.20*a[15][t]+26923086.46*a[17][t]
37522329.49*a[17][t]-.1082338049e11*a[18][t]-27435155.86*a[16][t]+9251036.956*a[13][t]-19358824.20*a[14][t]
-.5529464441e11*a[18][t]-27435155.86*a[15][t]-20873274.82*a[14][t]-6275707.162*a[13][t]-.1158227844e11*a[17][t]
-.1615478137e12*a[18][t]+26923086.46*a[14][t]+37522329.49*a[15][t]-5539558.102*a[13][t]-.1158227844e11*a[16][t]
-.1082338049e11*a[15][t]-.5529464441e11*a[16][t]-.1615478137e12*a[17][t]+23266721.99*a[13][t]+25978444.80*a[14][t]

as u see,six terms are not included in the matrix A :
.1689228391e12*a[13][t]^2+6406045.412*a[14][t]^2-4317791.317*a[15][t]^2+2291210.983*a[16][t]^2-.2842690977e11*a[17][t]^2-.1782456232e12*a[18][t]^2

what should i do? 

could u  rewrite the code after edition?i could not understand what u mean by this edit?,
this code has another problem also, in matrix A  some expressions are not included,the expression including the power of 2 of the function,for example :
eq1:=.1689228391e12*a[13][t]^2+6406045.412*a[14][t]^2-4317791.317*a[15][t]^2+2291210.983*a[16][t]^2-.2842690977e11*a[17][t]^2-.1782456232e12*a[18][t]^2+.9846526429e12*a[1][t]+.2533881291e12*a[2][t]+.3076607771e11*a[3][t]+8105875.203*a[4][t]-5054889.363*a[5][t]-34561.30764*a[6][t]+9251036.956*a[15][t]*a[13][t]-19358824.20*a[14][t]*a[15][t]+.4614468816e11*a[13][t]*a[14][t]-.1082338049e11*a[18][t]*a[15][t]+23266721.99*a[18][t]*a[13][t]+25978444.80*a[18][t]*a[14][t]-.5529464441e11*a[18][t]*a[16][t]-.1615478137e12*a[17][t]*a[18][t]-27435155.86*a[16][t]*a[15][t]-6275707.162*a[16][t]*a[13][t]-20873274.82*a[14][t]*a[16][t]+37522329.49*a[15][t]*a[17][t]-5539558.102*a[17][t]*a[13][t]+26923086.46*a[14][t]*a[17][t]-.1158227844e11*a[16][t]*a[17][t]:

 but the terms of A(1,i) i=1..18 are :

.9846526429e12
.2533881291e12
.3076607771e11
0.8105875203e7
-.5054889363e7
-34561.30764
0
0
0
0
0
0
23266721.99*a[18][t]-6275707.162*a[16][t]+.4614468816e11*a[14][t]-5539558.102*a[17][t]+9251036.956*a[15][t]
.4614468816e11*a[13][t]-20873274.82*a[16][t]+25978444.80*a[18][t]-19358824.20*a[15][t]+26923086.46*a[17][t]
37522329.49*a[17][t]-.1082338049e11*a[18][t]-27435155.86*a[16][t]+9251036.956*a[13][t]-19358824.20*a[14][t]
-.5529464441e11*a[18][t]-27435155.86*a[15][t]-20873274.82*a[14][t]-6275707.162*a[13][t]-.1158227844e11*a[17][t]
-.1615478137e12*a[18][t]+26923086.46*a[14][t]+37522329.49*a[15][t]-5539558.102*a[13][t]-.1158227844e11*a[16][t]
-.1082338049e11*a[15][t]-.5529464441e11*a[16][t]-.1615478137e12*a[17][t]+23266721.99*a[13][t]+25978444.80*a[14][t]

as u see,six terms are not included in the matrix A :
.1689228391e12*a[13][t]^2+6406045.412*a[14][t]^2-4317791.317*a[15][t]^2+2291210.983*a[16][t]^2-.2842690977e11*a[17][t]^2-.1782456232e12*a[18][t]^2

what should i do? 

@Carl Love 
i have a problem,

M:= dsolve(SYS union inits,numeric, method = dverk78,tolerance=1);
> M(0.000001);
M := proc(x_dverk78) ... end proc

[t = 0.100000000000000 10e-5 , a[1, 1](t) = -0.959988511261622 e10 ,

d
-- a[1, 1](t) = -13.5037741005812,
dt

a[1, 2](t) = 2.91235173166318,

d 7
-- a[1, 2](t) = 0.936180472200254 e10 ,
dt

a[1, 3](t) = 0.0144117978131728,

d
-- a[1, 3](t) = 14801.1119104049,
dt

a[2, 1](t) = 0.000240937360366017,

d
-- a[2, 1](t) = 417.905960165959,
dt

a[2, 2](t) = -15.3372441258150,

d 8
-- a[2, 2](t) = -0.496862993966895 e10 ,
dt

a[2, 3](t) = -0.0229010942815061,

d
-- a[2, 3](t) = -112758.921973661,
dt

a[3, 1](t) = -0.00128534036908000,

d
-- a[3, 1](t) = -2256.83712143968,
dt

a[3, 2](t) = 45.8168205401814,

d 9
-- a[3, 2](t) = 0.148794637685973 e10 ,
dt

a[3, 3](t) = 0.116389577971581,

d
-- a[3, 3](t) = 371189.882914135,
dt

a[4, 1](t) = 0.00386564225786822,

d
-- a[4, 1](t) = 7078.52931205466,
dt

a[4, 2](t) = -67.9420540103526,

d 9
-- a[4, 2](t) = -0.221188329407724 e10 ,
dt

a[4, 3](t) = -0.157827175344053,

d
-- a[4, 3](t) = -597031.914761630,
dt

a[5, 1](t) = -0.00401999280068150,

d
-- a[5, 1](t) = -7410.47508303779,
dt

a[5, 2](t) = 52.0149581038077,

d 9
-- a[5, 2](t) = 0.169615117314704 e10 ,
dt

a[5, 3](t) = 0.139243986965477,

d
-- a[5, 3](t) = 469925.748346270,
dt

a[6, 1](t) = 0.0121358461805532,

d
-- a[6, 1](t) = 4057.75033209267,
dt

a[6, 2](t) = -15.2004439836840,

d 8
-- a[6, 2](t) = -0.496906458911828 e10 ,
dt

a[6, 3](t) = -0.0280377065461516,

d
-- a[6, 3](t) = -142529.037267704]
dt
when i change the options of soving,i face this :

M:= dsolve(SYS union inits,numeric, method = dverk78, abserr = 1.*10^(-8), relerr = 1.*10^(-8));
> M(0.000001);

M := proc(x_dverk78) ... end proc
Error, (in M) cannot evaluate the solution past 0., step size < hmin, problem may be singular or error tolerance may be too small

 i change the abserr to 1 :

M:= dsolve(SYS union inits,numeric, method = dverk78, abserr = 1, relerr =1);
> M(0.000001);
M := proc(x_dverk78) ... end proc

-5 -5
[t = 0.100000000000000 10 , a[1, 1](t) = -0.959988511261595 10 ,

d
-- a[1, 1](t) = -13.5037741005798,
dt

a[1, 2](t) = 2.91235173166317,

d 7
-- a[1, 2](t) = 0.936180472200238 10 ,
dt

a[1, 3](t) = 0.0144117978131727,

d
-- a[1, 3](t) = 14801.1119104041,
dt

a[2, 1](t) = 0.000240937360366015,

d
-- a[2, 1](t) = 417.905960165948,
dt

a[2, 2](t) = -15.3372441258150,

d 8
-- a[2, 2](t) = -0.496862993966902 10 ,
dt

a[2, 3](t) = -0.0229010942815055,

d
-- a[2, 3](t) = -112758.921973657,
dt

a[3, 1](t) = -0.00128534036907998,

d
-- a[3, 1](t) = -2256.83712143960,
dt

a[3, 2](t) = 45.8168205401811,

d 9
-- a[3, 2](t) = 0.148794637685973 10 ,
dt

a[3, 3](t) = 0.116389577971577,

d
-- a[3, 3](t) = 371189.882914122,
dt

a[4, 1](t) = 0.00386564225786819,

d
-- a[4, 1](t) = 7078.52931205448,
dt

a[4, 2](t) = -67.9420540103533,

d 9
-- a[4, 2](t) = -0.221188329407726 10 ,
dt

a[4, 3](t) = -0.157827175344047,

d
-- a[4, 3](t) = -597031.914761610,
dt

a[5, 1](t) = -0.00401999280068146,

d
-- a[5, 1](t) = -7410.47508303753,
dt

a[5, 2](t) = 52.0149581038074,

d 9
-- a[5, 2](t) = 0.169615117314706 10 ,
dt

a[5, 3](t) = 0.139243986965471,

d
-- a[5, 3](t) = 469925.748346254,
dt

a[6, 1](t) = 0.0121358461805532,

d
-- a[6, 1](t) = 4057.75033209254,
dt

a[6, 2](t) = -15.2004439836838,

d 8
-- a[6, 2](t) = -0.496906458911823 10 ,
dt

a[6, 3](t) = -0.0280377065461499,

d
-- a[6, 3](t) = -142529.037267700]
dt

but when i change the time from 1e-6 to 1e-5 i face this :

M := proc(x_dverk78) ... end proc

Error, (in M) powering may produce overflow

 this means i have a singularity between the time 1-e5 and 1e-6,what should i do?can i do anything to my equations to avoid these singularites?

when i change the range of solving,i face this problem:

M:= dsolve(SYS union inits,numeric, stiff=true,range=1..10);Warning, cannot evaluate the solution further right of .43025197e-5, maxfun limit exceeded (see ?dsolve,maxfun for details)

                M := proc(x_rosenbrock)  ...  end proc 

M(2);
Error, (in M) cannot evaluate the solution further right of .43025197e-5, maxfun limit exceeded (see ?dsolve,maxfun for details)
and thus it has been solved for the range=1..10,but i can not see the results :( .
what should i do?
thnx for your help. 

yes,i have tried it numerically, but it does not answered goog,my equations and the discussion related are in this post :
http://www.mapleprimes.com/questions/147322-Dsolvenumeric-Takes-More-Than-12-Hours
thnx,

yes,i have tried it numerically, but it does not answered goog,my equations and the discussion related are in this post :
http://www.mapleprimes.com/questions/147322-Dsolvenumeric-Takes-More-Than-12-Hours
thnx,

@Carl Love i have 18 unknown function in my system of odes,and i waant to use step by step method to solve it,but i can not define neither iniitial conditon nor a loop for it?

i have complicated set of odes,and dsolve does not help me to solve them,thus i want to use step by step method to solve my odes, but i can not define a loop for my system of odes to solve , how can i define a loop for a system of ode?