Question: How should I present my system of first order ODE's to Maple to find a solution?

Hello,

I have a question regarding the dsolve function in Maple. I am trying to solve a system of 5 first order ODE's. First I solve the differential equations to arrive at the general solution using the dsolve command. Here Maple already produces very large and slow output while I would actually expect a much more compact output. 

Then I want to evaluate the final solution by using the initial conditons. At this point Maple keeps 'evaluating..' and does not come up with a solution. Only when I fill in all the parameters Maple finds a solution ( still very slow ). However as I want to be able to play with the input parameters after evaluation I want a solution without having to fill in the parameters prior to solving the system.

My question is whether I am presenting the system of ODE's in the right way to Maple. Should I rewrite the system for Maple to be able to produce a more dense solution? Or should I for instance use Laplace Transform to simplify the equations prior to solving?

Please help!

Joa

restart;

eq1 := sig_total-sig2(t)-sig3(t)-sig4(t)-sig5(t)+T1*(-(diff(sig2(t), t))-(diff(sig3(t), t))-(diff(sig4(t), t))-(diff(sig5(t), t))) = u1*(diff(eps(t), t));
eq2 := sig2(t)+T2*(diff(sig2(t), t)) = u2*(diff(eps(t), t));
eq3 := sig3(t)+T3*(diff(sig3(t), t)) = u3*(diff(eps(t), t));
eq4 := sig4(t)+T4*(diff(sig4(t), t)) = u4*(diff(eps(t), t));
eq5 := sig5(t)+T5*(diff(sig5(t), t)) = u5*(diff(eps(t), t));

dsolve({eq1, eq2, eq3, eq4, eq5}, {eps(t), sig2(t), sig3(t), sig4(t), sig5(t)});
sig_total := sig_0;
desys := {eq1, eq2, eq3, eq4, eq5}; ic := {eps(0) = sig_0/(E1+E2+E3+E4+E5), sig2(0) = sig_0/(1+(E1+E3+E4+E5)/E2), sig3(0) = sig_0/(1+(E1+E2+E4+E5)/E3), sig4(0) = sig_0/(1+(E1+E3+E2+E5)/E4), sig5(0) = sig_0/(1+(E1+E3+E4+E2)/E5)};

solution := combine(dsolve(desys union ic, {eps(t), sig2(t), sig3(t), sig4(t), sig5(t)})); assign(solution);

E_total := 33112; a := .1; b := .2; c := .15; d := .2; e := .35;

E1 := a*E_total; E2 := b*E_total; E3 := c*E_total; E4 := d*E_total; E5 := e*E_total; T1 := 1; T2 := 10; T3 := 100; T4 := 1000; T5 := 10000; u1 := T1*E1; u2 := T2*E2; u3 := T3*E3; u4 := T4*E4; u5 := T5*E5; sig_0 := 2;

plot(eps(t), t = 0 .. 672, y = 0 .. 0.12e-3);
y := eval(eps(t), t = 672);
evalf(y);
 

Download Maxwell_solve_new_method_5_units.mwMaxwell_solve_new_method_5_units.mw