Question: dsolve taking forever to solve?

Is there any way to simplify this code?  even with small number at discretization this take forever to solve  :(

restart;

with(plots);

numero := 5;

# Valores Calculados / Retirados da internet

viscosidade := Vector[row](10, [0.000024, 0.00001113, 0.89*10^(-5), 0.00001779, 0.017299, 0.00001028, 0.927*10^(-5), 0.818*10^(-5), 0.749*10^(-5), 0.7*10^(-5)]);

H := Vector[row](8, [-0.1414243148*10^8, -0.1677843875*10^8, -0.2177577229*10^8, -0.3078557121*10^8, -0.4007777822*10^8, -0.4007777822*10^8, -0.5832487417*10^8, 0.84247483*10^7]);

a := [524, 879, 1271, 1099, 1779, -163, 241, -258, 1200, 1200];

b := [1.3383, 4.1117, 0.8467, -0.47, 0.0333, 6.78, 5.6683, 7.5933, 3, 3];
c := [-0.0008, 0.00015, -0.00145, 0.00105, 0.0008, -0.00475, -0.002, -0.00455, 0, 0];
d := [0.166667*10^(-6), -0.666667*10^(-6), 0.833333*10^(-6), -0.5*10^(-6), -0.333333*10^(-6), 0.15*10^(-5), 0.166667*10^(-6), 0.116667*10^(-5), -0.888178*10^(-15), -0.888178*10^(-15)];

Massas := Vector[row](10, [44.01, 16.04, 2.02, 28.01, 18.01, 28.05, 30.07, 44.1, 58.12, 84]);

nu := Matrix(11, 8, [[-1, -2, -2, -3, -4, -4, -6.05, -1], [-3, -4, -5, -7, -9, -9, -12.23, 1], [1, 2, 2, 3, 4, 4, 6.05, -1], [1, 0, 0, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0, 0, 0], [0, 0, 1, 0, 0, 0, 0, 0], [0, 0, 0, 1, 0, 0, 0, 0], [0, 0, 0, 0, 1, 0, 0, 0], [0, 0, 0, 0, 0, 1, 0, 0], [0, 0, 0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 0, 0, 0, 1]]);

A := Vector[row](8, [65.56165, 0.045, 0.001144, 0.119*10^(-5), 0.184*10^(-9), 0.659*10^(-8), 0.198*10^(-7), 0.0000378]);
E := Vector[row](8, [83523.900, 65017, 49782, 34885.500, 27728.900, 25730.100, 23564.300, 58826.300]);
m := Vector[row](8, [0.60950, 0.37330, -0.02830, 0.43540, 0.03530, -0.25150, 0.81870, -0.36200]);
n := Vector[row](8, [-0.46790, -0.26130, 0.19130, 0.08760, 1.11630, 0.00920, 0.03420, 1.26390]);

Conc[CO] := 40.65;
Conc[H2] := 78.91;
Conc[N2] := 119.5;

DensidadeCO := 6.712;
DensidadeN2 := 6.736;
DensidadeH2 := 0.0483;

Cpcte := [1.064, 14.51, 1.954, 2.889, 2.227, 2.596, 2.583, 2.619, 2.619, 2.619, 0.856, 1.056];
viscosidade := Vector[row](12, [0.0000252, 0.0000252, 0.0000252, 0.0000163, 0.0000154, 0.0000154, 0.0000154, 0.0000154, 0.0000154, 0.000023, 0.000023, 0.0000251]);
massas := Vector[row](12, [0.02801, 0.00202, 0.01801, 0.01604, 0.02805, 0.03007, 0.0441, 0.05812, 0.05812, 0.084, 0.04401, 0.02801]);

#Constantes
Rcte := 8.314;
varepsilon[B] := 0.4;
beta := 0.255;
phi := 4/3;
delta[wall] := 0.004;
lambda[wall] := 60;
h[water] := 1000;
Temp[water] := 498;
rho[B] := 380;
dt := 0.0157;
dp := 0.00015;

#Equações Auxiliares
for k to 11 do
    Y[k](t, z) := Rcte*C[k](t, z)*T(t, z)/1000000;
end do;

Cp[mix] := (t, z) -> 1000*Cpcte[12];
mu[mix] := (t, z) -> viscosidade[12];
M[mix] := (t, z) -> massas[12];

h[int] := (t, z) -> rho[mix](t, z)*0.458/varepsilon[B]/((Cp[mix](t, z)*mu[mix](t, z))^4.074*(dp/mu[mix](t, z))^4.407);
Reynolds := (t, z) -> dp*rho[mix](t, z)/mu[mix](t, z);
f := (t, z) -> 172/Reynolds(t, z) + 4.36/Reynolds(t, z)^0.12;
U := (t, z) -> 1/(1/h[water] + delta[wall]/lambda[wall] + 1/h[int](t, z));
P[CO] := (t, z) -> 1000000*Y[1](t, z);
P[H2] := (t, z) -> 1000000*Y[2](t, z);
for j to 8 do
    R[j](t, z) := A[j]*exp(-E[j]/(Rcte*T(t, z)))*P[CO](t, z)^m[j]*P[H2](t, z)^n[j];
end do;

#EDP's

for k to 11 do
    edp[k] := diff(C[k](t, z), t) = -v(t, z)*diff(C[k](t, z), z) + rho[B]*beta*add(nu[k][j]*R[j](t, z), j = 1 .. 8);
end do;
edp[12] := diff(T(t, z), t) = -v(t, z)*diff(T(t, z), z) + rho[B]*beta*add(add(-H[j]*nu[i][j]*R[j](t, z), i = 1 .. 10), j = 1 .. 8)/(rho[mix](t, z)*Cp[mix](t, z)) - 4*U(t, z)*(T(t, z) - Temp[water])/(dt*rho[mix](t, z)*Cp[mix](t, z));
edp[13] := diff(v(t, z), z) = -v(t, z)*diff(rho[mix](t, z), z)/rho[mix](t, z);
edp[14] := diff(rho[mix](t, z), z) = M[mix](t, z)/Rcte*(-1000000*diff(T(t, z), z)/T(t, z)^2);
edp[15] := diff(PT(t, z), z) = -f(t, z)*v(t, z)^2*rho[mix](t, z)/dp;

#Discretização do Modelo
hh := 0.11/numero;
for k to 11 do
    dis[k] := diff(C[k](t, z), z) = (x[k, i](t) - x[k, i - 1](t))/hh;
end do;
dis[12] := diff(T(t, z), z) = (x[12, i](t) - x[12, i - 1](t))/hh;
dis[13] := diff(rho[mix](t, z), z) = (x[13, i](t) - x[13, i - 1](t))/hh;
dis[14] := diff(v(t, z), z) = (x[14, i](t) - x[14, i - 1](t))/hh;
dis[15] := diff(PT(t, z), z) = (x[15, i](t) - x[15, i - 1](t))/hh;
for k from 16 to 26 do
    dis[k] := C[k - 15](t, z) = x[k - 15, i](t);
end do;
dis[27] := T(t, z) = x[12, i](t);
dis[28] := rho[mix](t, z) = x[13, i](t);
dis[29] := v(t, z) = x[14, i](t);
dis[30] := PT(t, z) = x[15, i](t);
listadis := seq(dis[k], k = 1 .. 30);
for k to 15 do
    equacao[k] := eval(edp[k], {listadis});
end do;

#Resolução 

ci[1, i] := x[1, i](0) = Conc[CO];
ci[2, i] := x[2, i](0) = Conc[H2];
for k from 3 to 11 do
    ci[k, i] := x[k, i](0) = 0;
end do;
ci[12, i] := x[12, i](0) = 503;
ci[13, i] := x[13, i](0) = 0.005165;
ci[14, i] := x[14, i](0) = 0.33*DensidadeH2 + 0.17*DensidadeCO + 0.5*DensidadeN2;
ci[15, i] := x[15, i](0) = 1000000;
cis := seq(seq(eval(ci[k, i], i = j), k = 1 .. 15), j = 1 .. numero);
unassign(i);
eqs := seq(seq(eval(equacao[k], i = j), k = 1 .. 15), j = 2 .. numero);
final[1] := x[1, 1](t) = Conc[CO];
final[2] := x[2, 1](t) = Conc[H2];
for k from 3 to 11 do
    final[k] := x[k, 1](t) = 0;
end do;
final[12] := x[12, 1](t) = 503;
final[13] := x[13, 1](t) = 0.005165;
final[14] := x[14, 1](t) = 0.33*DensidadeH2 + 0.17*DensidadeCO + 0.5*DensidadeN2;
final[15] := x[15, 1](t) = 1000000;
seqfinal := seq(final[k], k = 1 .. 15), eqs;
sol := dsolve({cis, seqfinal}, numeric, stiff = true, range = 0 .. 180);
odeplot(sol, [t, x[12, 4](t)], t = 0 .. 180);

Modelo_discretizado_com_t_e_tudo_constante.mw

Thank you already