solve a system of 8 coupled ODEs

Question

I'm trying to solve this differential equation system in matlab:

with this code:

syms w1(t) w1c(t) w2(t) w2c(t) m1(t) m1c(t) m2(t) m2c(t) d a1w b1w y1w a2w b2w y2w a1m b1m y1m a2m b2m y2m
ode1 = diff(w1,t) == (-d*w1)-(a1w*w1)-(b1w*(w1*(w1c+m1c)))+(y1w*w1);
ode2 = diff(w1c,t) == (a1w*w1)+(b1w*(w1*(w1c+m1c)))-(y1w*w1c);
ode3 = diff(w2,t) == (d*w1)-(a2w*w2)-(b2w*(w2*(w2c+m2c)))+(y2w*w2c);
ode4 = diff(w2c,t) == (a2w*w2)+(b2w*(w2*(w2c+m2c)))-(y2w*w2c);
ode5 = diff(m1,t)==-(d*m1)-(a1m*m1)-(b1m*(m1*(w1c+m1c)))+(y1m*m1c);
ode6 = diff(m1c,t) == (a1m*m1)+(b1m*(m1*(w1c+m1c)))-(y1m*m1c);
ode7 = diff(m2,t) == (d*m1)-(a2m*m2)-(b2m*(m2*(w2c+m2c)))+(y2m*m2c);
ode8 = diff(m2c,t) == (a2m*m2)+(b2m*(m2*(w2c+m2c)))-(y2m*m2c);

odes = [ode1; ode2; ode3; ode4; ode5; ode6; ode7; ode8]
S = dsolve(odes)

with this initial conditions:

cond1 = w1(0) == 19.86;
cond2 = w2(0) == 16.21;
cond3 = w1c(0) == 9.75;
cond4 = w2c(0) == 9.52;
cond5 = m1(0) == 15.24;
cond6 = m2(0) == 11.64;
cond7 = m1c(0) == 9.33;
cond8 = m2c(0) == 8.63;
conds = [cond1; cond2; cond3; cond4; cond5; cond6; cond7; cond8];
[w1Sol(t), w1cSol(t), w2Sol(t), w2cSol(t), m1Sol(t), m1cSol(t), m2Sol(t), m2cSol(t)] = dsolve(odes,conds)

but i get this text:

Warning: Explicit solution could not be found. 
> In dsolve at 197 

S =

[ empty sym ]

I've read about the ode45 function but I don't know how to use it properly so this is what I did:

f=@(t,x) [(-d*x(1))-(a1w*x(1))-(b1w*(x(1)*(x(2)+x(6))))+(y1w*x(1));
    (a1w*x(1))+(b1w*(x(1)*(x(2)+x(6))))-(y1w*x(2));
    (d*x(1))-(a2w*x(3))-(b2w*(x(3)*(x(4)+x(8))))+(y2w*x(4));
    (a2w*x(3))+(b2w*(x(3)*(x(4)+x(8))))-(y2w*x(4));
    -(d*x(5))-(a1m*x(5))-(b1m*(x(5)*(x(2)+x(6))))+(y1m*x(6));
    (a1m*x(5))+(b1m*(x(5)*(x(2)+x(6))))-(y1m*x(6));
    (d*x(5))-(a2m*x(7))-(b2m*(x(7)*(x(4)+x(8))))+(y2m*x(8));
    (a2m*x(7))+(b2m*(x(7)*(x(4)+x(8))))-(y2m*x(8))]

I don't know if it's right or now, or how to put my initial conditions there. Any help will be appreciated.

UPDATE: I put the defined parameters as

a1w = 0.0180156
a1m = 0.0180175
a2w = 0.000119175
a2m = 0.00307025
b1w = 0.538867
b1m = 2.59115
b2w = 0.668998
b2m = 1.38138
y1w = 0.285557
y1m = 0.853703
y2w = 0.270119
y2m = 0.405319

and run with this:

[t,xa] = ode45(f,[0 8],[19.86 16.21 9.75 9.52 15.24 11.64 9.33 8.63]);

Am I doing this correctly? After this I'm supposed to plot every one but also i'm not very sure how.