I want to solve the following Matrix equation:

![enter image description here][1]


  [1]: https://i.sstatic.net/h7ZaL.png

I have deffined everything in the following code:

    r24 =.;
    r24[s24_] = {{1, 0, 0, 0}, {0, Sqrt[1 - s24^2], 0, s24}, {0, 0, 1,0}, {0,-s24, 0, Sqrt[1- s24^2]}};


    r23 = {{1, 0, 0, 0}, {0, 0.707, 0.707, 0}, {0, -0.707, 0.707, 0}, {0, 0, 0, 1}};
    r13 = {{0.948, 0, 0.32, 0}, {0, 1, 0, 0}, {-0.32, 0, 0.948, 0}, {0, 0,0, 1}};
    r12 = {{0.837, 0.548, 0, 0}, {-0.548, 0.837, 0, 0}, {0, 0, 1, 0}, {0, 0, 0, 1}};

    u =.;
    u[s24_] := r24[s24].r23.r13.r12;

    udager =.;
    udager[s24_] := Transpose[u[s24]];

    m2 =.;
    m2[w41_] := {{0, 0, 0, 0}, {0, 8*10^-5, 0, 0}, {0, 0, 2.32*10^-3,0}, {0, 0, 0, w41}};

    term1 =.;
    term1[F_, s24_, w41_] := (1/2 F)*(u[s24].m2[w41].udager[s24]);

    term1[F, s24, w41] // MatrixForm;

    a = 6371;

    r[x_, cosz_] := Sqrt[x^2 + a^2 - 2*a*x*cosz]

    ne[x_, cosz_] :=Piecewise[{{13.088 - 8.84*(r[x, cosz]/a)^2,0 <= r[x, cosz] < 1221.5}, {12.58 - 1.26*(r[x, cosz]/a) -3.64*(r[x, cosz]/a)^2 - 5.53*(r[x, cosz]/a)^3, 
    1221.5 <= r[x, cosz] < 3480}, {7.96 - 6.48*(r[x, cosz]/a) + 5.53*(r[x, cosz]/a)^2 - 3.08*(r[x, cosz]/a)^3,3480 <= r[x, cosz] < 3630}, {7.96 - 6.48*(r[x, cosz]/a) + 
    5.53*(r[x, cosz]/a)^2 - 3.08*(r[x, cosz]/a)^3,3630 <= r[x, cosz] < 5600}, {7.96 - 6.48*(r[x, cosz]/a) +5.53*(r[x, cosz]/a)^2 - 3.08*(r[x, cosz]/a)^3, 
    5600 <= r[x, cosz] < 5701}, {5.32 - 1.48*(r[x, cosz]/a), 
    5701 <= r[x, cosz] < 5771}, {11.25 - 8.03*(r[x, cosz]/a), 
    5771 <= r[x, cosz] < 5971}, {7.11 - 3.8*(r[x, cosz]/a), 
    5971 <= r[x, cosz] < 6151}, {2.69 + 0.69*(r[x, cosz]/a), 
    6151 <= r[x, cosz] < 6291}, {2.69 + 0.69*(r[x, cosz]/a), 
    6291 <= r[x, cosz] < 6346.6}, {2.9,6346.6 <= r[x, cosz] < 6356}, {2.6,6356 <= r[x, cosz] < 6368}, {1.02, 6368 <= r[x, cosz] < 6371}}]

    vr[x_, cosz_] :=7.63*10^(-14)*{{ne[x, cosz], 0, 0, 0}, {0, 0, 0, 0}, {0, 0, 0,0}, {0, 0, 0, ne[x, cosz]/2}};

    psi[x_, cosz_, w41_, s24_, F_] := {psi1[x, cosz, w41, s24, F],psi2[x, cosz, w41, s24, F], psi3[x, cosz, w41, s24, F],psi4[x, cosz, w41, s24, F]}

    Derivative[{1, 0, 0, 0, 0}][psi][{x, cosz, w41, s24,F}] = (1/Sqrt[-1])*((term1[F, s24, w41] + vr[x, cosz]).psi[x, cosz, w41,s24, F]);

With the following initial conditions:  

    psi1[-2*a*cosz, cosz, w41, s24, F] = 0; psi2[-2*a*cosz, cosz, w41, s24, F] = 1; 
    psi3[-2*a*cosz, cosz, w41, s24, F] = 0; psi4[-2*a*cosz, cosz, w41, s24, F] = 0;

and `x` changes beetween `-2*a*cosz` and `0`. My goal is to find `psi2[cosz, w41, s24,F]`. 

I can't just use `DSolve`, because it is not possible to solve such this system of equations analytically. And I don't know how I can use something like `NDSolve` when I have these many parameters. Can any body help me with this?