# System of 16 coupled Integro-differential equations

I have a question: How can i solve numerically a system of 16 coupled integro-differential equations?

$$- a{{s'}_{11}} + b{{s'}_{33}} + c{{s'}_{44}} + Ng_1^2(\int_0^z {({s_{14}}(t) + {s_{23}}(t))} dt){{s'}_{41}} + Ng_1^2(\int_0^z {({s_{41}}(t) + {s_{32}}(t))} dt){{s'}_{14}} + Ng_2^2(\int_0^z {({s_{13}}(t) + {s_{24}}(t))} dt){{s'}_{31}} + Ng_2^2(\int_0^z {({s_{31}}(t) + {s_{42}}(t))} dt){{s'}_{13}} = - (Ng_1^2 + Ng_2^2){s_{11}} - Ng_2^2{s_{33}} - Ng_1^2{s_{44}}\\$$ $$- a{{s'}_{22}} + b{{s'}_{44}} + c{{s'}_{33}} + Ng_1^2(\int_0^z {({s_{14}}(t) + {s_{23}}(t))} dt){{s'}_{32}} + Ng_1^2(\int_0^z {({s_{41}}(t) + {s_{32}}(t))} dt){{s'}_{23}} + Ng_2^2(\int_0^z {({s_{13}}(t) + {s_{24}}(t))} dt){{s'}_{42}} + Ng_2^2(\int_0^z {({s_{31}}(t) + {s_{42}}(t))} dt){{s'}_{24}} = - (Ng_1^2 + Ng_2^2){s_{22}} - Ng_2^2{s_{44}} - Ng_1^2{s_{33}}\\$$ $$- (a + d){{s'}_{33}} - Ng_1^2(\int_0^z {({s_{14}}(t) + {s_{23}}(t))} dt){{s'}_{32}} - Ng_1^2(\int_0^z {({s_{41}}(t) + {s_{32}}(t))} dt){{s'}_{23}} - Ng_2^2(\int_0^z {({s_{13}}(t) + {s_{24}}(t))} dt){{s'}_{31}} - Ng_2^2(\int_0^z {({s_{31}}(t) + {s_{42}}(t))} dt){{s'}_{13}} = - (Ng_1^2 + Ng_2^2){s_{33}} + Ng_2^2{s_{11}} + Ng_1^2{s_{22}}\\$$ $$- (a + d){{s'}_{44}} - Ng_1^2(\int_0^z {({s_{14}}(t) + {s_{23}}(t))} dt){{s'}_{41}} - Ng_1^2(\int_0^z {({s_{41}}(t) + {s_{32}}(t))} dt){{s'}_{14}} - Ng_2^2(\int_0^z {({s_{13}}(t) + {s_{24}}(t))} dt){{s'}_{42}} - Ng_2^2(\int_0^z {({s_{31}}(t) + {s_{42}}(t))} dt){{s'}_{24}} = - (Ng_1^2 + Ng_2^2){s_{44}} + Ng_2^2{s_{22}} + Ng_1^2{s_{11}}\\$$ $${a_1}{{s'}_{41}} - Ng_1^2(\int_0^z {({s_{41}}(t) + {s_{32}}(t))} dt)({{s'}_{11}} - {{s'}_{44}}) - Ng_2^2(\int_0^z {({s_{31}}(t) + {s_{42}}(t))} dt)({{s'}_{21}} - {{s'}_{43}}) = (Ng_1^2 + Ng_2^2){s_{41}}\\$$ $${a_2}{{s'}_{14}} + Ng_1^2(\int_0^z {({s_{14}}(t) + {s_{23}}(t))} dt)({{s'}_{44}} - {{s'}_{11}}) + Ng_2^2(\int_0^z {({s_{13}}(t) + {s_{24}}(t))} dt)({{s'}_{34}} - {{s'}_{21}}) = - (Ng_1^2 + Ng_2^2){s_{14}}\\$$ $${b_1}{{s'}_{31}} - Ng_1^2(\int_0^z {({s_{41}}(t) + {s_{32}}(t))} dt)({{s'}_{21}} - {{s'}_{34}}) - Ng_2^2(\int_0^z {({s_{31}}(t) + {s_{42}}(t))} dt)({{s'}_{11}} - {{s'}_{33}}) = (Ng_1^2 + Ng_2^2){s_{31}}\\$$ $${b_2}{{s'}_{14}} + Ng_1^2(\int_0^z {({s_{14}}(t) + {s_{23}}(t))} dt)({{s'}_{43}} - {{s'}_{12}}) + Ng_2^2(\int_0^z {({s_{13}}(t) + {s_{24}}(t))} dt)({{s'}_{33}} - {{s'}_{11}}) = - (Ng_1^2 + Ng_2^2){s_{13}}\\$$ $${b_1}{{s'}_{42}} - Ng_1^2(\int_0^z {({s_{41}}(t) + {s_{32}}(t))} dt)({{s'}_{12}} - {{s'}_{43}}) - Ng_2^2(\int_0^z {({s_{31}}(t) + {s_{42}}(t))} dt)({{s'}_{22}} - {{s'}_{44}}) = (Ng_1^2 + Ng_2^2){s_{42}}\\$$ $${b_2}{{s'}_{24}} + Ng_1^2(\int_0^z {({s_{14}}(t) + {s_{23}}(t))} dt)({{s'}_{34}} - {{s'}_{12}}) + Ng_2^2(\int_0^z {({s_{13}}(t) + {s_{24}}(t))} dt)({{s'}_{44}} - {{s'}_{22}}) = - (Ng_1^2 + Ng_2^2){s_{24}}\\$$ $${c_1}{{s'}_{32}} - Ng_1^2(\int_0^z {({s_{41}}(t) + {s_{32}}(t))} dt)({{s'}_{22}} - {{s'}_{33}}) - Ng_2^2(\int_0^z {({s_{31}}(t) + {s_{42}}(t))} dt)({{s'}_{21}} - {{s'}_{34}}) = (Ng_1^2 + Ng_2^2){s_{32}}\\$$ $${c_2}{{s'}_{24}} + Ng_1^2(\int_0^z {({s_{14}}(t) + {s_{23}}(t))} dt)({{s'}_{33}} - {{s'}_{22}}) + Ng_2^2(\int_0^z {({s_{13}}(t) + {s_{24}}(t))} dt)({{s'}_{43}} - {{s'}_{21}}) = - (Ng_1^2 + Ng_2^2){s_{23}}\\$$ $$- (a + d){{s'}_{43}} - Ng_1^2(\int_0^z {({s_{14}}(t) + {s_{23}}(t))} dt){{s'}_{42}} - Ng_1^2(\int_0^z {({s_{41}}(t) + {s_{32}}(t))} dt){{s'}_{13}} - Ng_2^2(\int_0^z {({s_{13}}(t) + {s_{24}}(t))} dt){{s'}_{41}} - Ng_2^2(\int_0^z {({s_{31}}(t) + {s_{42}}(t))} dt){{s'}_{23}} = (Ng_1^2 + Ng_2^2){s_{43}} + Ng_2^2{s_{21}} + Ng_1^2{s_{12}}\\$$ $$- (a + d){{s'}_{34}} - Ng_1^2(\int_0^z {({s_{14}}(t) + {s_{23}}(t))} dt){{s'}_{31}} - Ng_1^2(\int_0^z {({s_{41}}(t) + {s_{32}}(t))} dt){{s'}_{24}} - Ng_2^2(\int_0^z {({s_{13}}(t) + {s_{24}}(t))} dt){{s'}_{32}} - Ng_2^2(\int_0^z {({s_{31}}(t) + {s_{42}}(t))} dt){{s'}_{14}} = (Ng_1^2 + Ng_2^2){s_{34}} + Ng_2^2{s_{12}} + Ng_1^2{s_{21}}\\$$ $$- a{{s'}_{21}} + Ng_1^2(\int_0^z {({s_{14}}(t) + {s_{23}}(t))} dt){{s'}_{31}} + Ng_1^2(\int_0^z {({s_{41}}(t) + {s_{32}}(t))} dt){{s'}_{24}} + Ng_2^2(\int_0^z {({s_{13}}(t) + {s_{24}}(t))} dt){{s'}_{41}} + Ng_2^2(\int_0^z {({s_{31}}(t) + {s_{42}}(t))} dt){{s'}_{23}} = - (Ng_1^2 + Ng_2^2){s_{21}} - (Ng_1^2 + Ng_2^2){s_{34}}\\$$ $$- a{{s'}_{12}} + Ng_1^2(\int_0^z {({s_{14}}(t) + {s_{23}}(t))} dt){{s'}_{42}} + Ng_1^2(\int_0^z {({s_{41}}(t) + {s_{32}}(t))} dt){{s'}_{31}} + Ng_2^2(\int_0^z {({s_{13}}(t) + {s_{24}}(t))} dt){{s'}_{32}} + Ng_2^2(\int_0^z {({s_{31}}(t) + {s_{42}}(t))} dt){{s'}_{14}} = - (Ng_1^2 + Ng_2^2){s_{12}} - Ng_1^2{s_{43}} + Ng_2^2{s_{34}}\\$$

we have 16 variable $${s_{ij}}$$ and they depend on z, we are looking for them.N, a,b,c,d,$${g_1}$$,$${g_2}$$,$${a_1}$$,$${a_2}$$,$${b_1}$$,$${b_2}$$,$${c_1}$$ and $${c_2}$$ are constant.

New contributor
Farhad.R is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
• Set $S_{ij}(z)=\int_0^zs_{ij}(s)\,ds$ so that your equations now are second order ODE in the $S_{ij}$. Do you see a point where that would not work? // It is often advisable to keep the original structure of the problem, it may also lead to a similarly structured implementation, helping to avoid "stupid" errors. – LutzL Aug 14 at 11:02
• So i don't get. whats your opinion to solve this problem? does Matlab or Mathematica help me or i have to write a code in C or python for example? – Farhad.R Aug 14 at 11:18
• I do not see clearly if the resulting system is then cleanly solvable for the second derivatives. If yes, it is an ODE system and can be solved by every ODE solver. Inform yourself on how the solver quality compares over the mentioned variants. If it is not solvable, it is a DAE system, might require index reduction, and is then solvable using DAE solvers. Many packages have this facility, however limited it may be. – LutzL Aug 14 at 11:30