I found a way to get a solution. By telling maple to use the Lapace method, it is able to find a symbolical solution to the problem. In the end my goal is to be able to fit the solution to the system of ODE's to selected data by varying the parameters a,b,c and d.

I am trying to fit the data using Directsearch, however I get an error:

My question thus is whether there is a way to fit data using the solution to the system of ODE's as presented in the Maple script. Please help!

Maple script:

restart;
Solving the system of differential equations with the initial conditions

eq1 := sig_total-sig2(t)-sig3(t)-sig4(t)+T1*(-(diff(sig2(t), t))-(diff(sig3(t), t))-(diff(sig4(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));

sig_total := sig_0;
desys := {eq1, eq2, eq3, eq4}; ic := {eps(0) = eps_0, sig2(0) = sig2_0, sig3(0) = sig3_0, sig4(0) = sig4_0};

solution := combine(dsolve(desys union ic, {eps(t), sig2(t), sig3(t), sig4(t)}, method = laplace)); assign(solution);
NULL;
Here you can play with the input parameters
E_total := 21943;
a := .25;
b := .25;
c := .2;
d := .3;

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

sig2_0 := sig_0/(1+(E1+E3+E4)/E2); sig3_0 := sig_0/(1+(E1+E2+E4)/E3); sig4_0 := sig_0/(1+(E1+E2+E3)/E4); eps_0 := sig_0/(E1+E2+E3+E4);

plot(eps(t), t = 0 .. 672, y = 0 .. 0.2e-3);

eps0 := sig_0/(E1+E2+E3+E4);
                    eps0 := 0.00004557261997
eps672 := eval(eps(t), t = 672);
evalf(eps672);
                        0.0001114003649

Here you can fit selected data X and Y (automatically?!)

a := 'a'; b := 'b'; c := 'c'; d := 'd';
E1 := a*E_total; E2 := b*E_total; E3 := c*E_total; E4 := d*E_total; T1 := 10; T2 := 100; T3 := 1000; T4 := 10000; u1 := T1*E1; u2 := T2*E2; u3 := T3*E3; u4 := T4*E4; sig_0 := 1;
sig2_0 := sig_0/(1+(E1+E3+E4)/E2); sig3_0 := sig_0/(1+(E1+E2+E4)/E3); sig4_0 := sig_0/(1+(E1+E2+E3)/E4); eps_0 := sig_0/(E1+E2+E3+E4);
NULL;
You can use the command above to show that eps(t) depends only on t and the parameters a,b,c and d.

epsfunc := eval(eps(t));

Provide the input data here:
X := Vector([0, 100, 200, 300, 400, 500], datatype = float);
Y := Vector([0.2e-2, 0.5e-2, 0.7e-2, 0.75e-2, 0.77e-2, 0.8e-2], datatype = float);
subplot1 := plot([0, 100, 200, 300, 400, 500], [0.2e-2, 0.5e-2, 0.7e-2, 0.75e-2, 0.77e-2, 0.8e-2]);

with(DirectSearch);
Digits := 15;
curvefit := DataFit(epsfunc, X, Y, t, [a = 0 .. 1, b = 0 .. 1, c = 0 .. 1, d = 0 .. 1, a+b+c+d = 1]);
F := eval(epsfunc, fit[2]);
subplot2 := plot(F, t = 0 .. 500, gridlines = true, linestyle = dashdot);
display({subplot1, subplot2}, size = [.4, .6]);
;

Maple_fitting_help.mw

@mmcdara Thank you for your help. However in your result the other condition which says that a,b and c should lie between 0..1 is no longer satisfied. Also, I have just edited the original post and added epsfunc.

@vv I edited the original post