% basic equation of material model WibCreep v5.1
function deps = odefun_materialmodel(t, eps, pE0, pE1, peta1, peta2, ...
    peta3, peta4, pE01d, peta11d, peta21d, peta31d, peta41d, sigma)

  % get model parameters
  E0    = e_aging(pE0, pE01d, t);
  E1    = zeros(length(t),1);
  E1(:) = pE1;
  E2 = E1;
  E3 = E1;
  E4 = E1;
  eta1 = eta_aging(peta1, peta11d, t);
  eta2 = eta_aging(peta2, peta21d, t);
  eta3 = eta_aging(peta3, peta31d, t);
  eta4 = eta_aging(peta4, peta41d, t);

  % estimate coefficients of ode
  a0 = 1/(E0*eta1*eta2*eta3*eta4);
  a1 = 1/(E0*E1*eta2*eta3*eta4) + 1/(E0*E2*eta1*eta3*eta4) + ...
      1/(E0*E3*eta1*eta2*eta4) + 1/(E0*E4*eta1*eta2*eta3);
  a2 = 1/(E0*E1*E2*eta3*eta4) + 1/(E0*E1*E3*eta2*eta4) + ...
      1/(E0*E1*E4*eta2*eta3) + 1/(E0*E2*E3*eta1*eta4) + ...
      1/(E0*E2*E4*eta1*eta3) + 1/(E0*E3*E4*eta1*eta2);
  a3 = 1/(E0*E1*E2*E3*eta4) + 1/(E0*E1*E2*E4*eta3) + ...
      1/(E0*E1*E3*E4*eta2) + 1/(E0*E2*E3*E4*eta1);
  a4 = 1/(E0*E1*E2*E3*E4);
  b0 = 1/(eta1*eta2*eta3*eta4); 
  b1 = 1/(E0*eta1*eta2*eta3) + 1/(E0*eta1*eta2*eta4) + ...
      1/(E0*eta1*eta3*eta4) + 1/(E0*eta2*eta3*eta4) + ...
      1/(E1*eta2*eta3*eta4) + 1/(E2*eta1*eta3*eta4) + ...
      1/(E3*eta1*eta2*eta4) + 1/(E4*eta1*eta2*eta3);
  b2 = 1/(E0*E1*eta2*eta3) + 1/(E0*E2*eta1*eta3) + ...
      1/(E0*E3*eta1*eta2) + 1/(E0*E1*eta2*eta4) + ...
      1/(E0*E2*eta1*eta4) + 1/(E0*E4*eta1*eta2) + 1/(E0*E1*eta3*eta4) + ...
      1/(E0*E3*eta1*eta4) + 1/(E0*E4*eta1*eta3) + 1/(E0*E2*eta3*eta4) + ...
      1/(E0*E3*eta2*eta4) + 1/(E0*E4*eta2*eta3) + 1/(E1*E2*eta3*eta4) + ...
      1/(E1*E3*eta2*eta4) + 1/(E1*E4*eta2*eta3) + 1/(E2*E3*eta1*eta4) + ...
      1/(E2*E4*eta1*eta3) + 1/(E3*E4*eta1*eta2);
  b3 = 1/(E0*E1*E2*eta3) + 1/(E0*E1*E3*eta2) + 1/(E0*E2*E3*eta1) + ...
      1/(E0*E1*E2*eta4) + 1/(E0*E1*E4*eta2) + 1/(E0*E2*E4*eta1) + ...
      1/(E0*E1*E3*eta4) + 1/(E0*E1*E4*eta3) + 1/(E0*E3*E4*eta1) + ...
      1/(E0*E2*E3*eta4) + 1/(E0*E2*E4*eta3) + 1/(E0*E3*E4*eta2) + ...
      1/(E1*E2*E3*eta4) + 1/(E1*E2*E4*eta3) + 1/(E1*E3*E4*eta2) + ...
      1/(E2*E3*E4*eta1);
  b4 = 1/(E0*E1*E2*E3) + 1/(E0*E1*E2*E4) + 1/(E0*E1*E3*E4) + ...
      1/(E0*E2*E3*E4) + 1/(E1*E2*E3*E4);

  % build system of ode
  deps = zeros(4,1);

  deps(1) = (a0*sigma-b3*eps(1)-b2*eps(2)-b1*eps(3)-b0*eps(4))/b4;
  deps(2) = eps(1);
  deps(3) = eps(2);
  deps(4) = eps(3);
end