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}}\\\tag1$$ $$ - 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}}\\\tag2$$ $$ - (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}}\\\tag3$$ $$ - (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}}\\\tag4$$ $${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}}\\\tag5$$ $${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}}\\\tag6$$ $${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}}\\\tag7$$ $${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}}\\\tag8$$ $${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}}\\\tag9$$ $${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}}\\\tag{10}$$ $${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}}\\\tag{11}$$ $${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}}\\\tag{12}$$ $$ - (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}}\\\tag{13}$$ $$ - (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}}\\\tag{14}$$ $$ - 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}}\\\tag{15}$$ $$ - 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}}\\\tag{16}$$
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.