Hello!

I am currently working on neutrino's oscillations and I need to solve a system of coupled ODE's. The problem is that I have never coded in Mathematica so my code is probably awful and... doesn't work at all. The system to solve is the following: (at the end, one must take the modulus of the amplitude)

with everything a constant except V_e a function of time. This is my code which does not return the right solution (see below)

Gf = 0;

Xmax = 100000;

O13 = N[Degree, 8.54];
O23 = N[Degree, 47.2];
O12 = N[Degree, 33.62];
dcp = N[Degree, 234];

dm21 = 7.40*10^(-5);
dm31 = 2.494*10^(-3);

p = 1;

M = {{0, 0, 0}, {0, dm21, 0}, {0, 0, dm31}};

R1 = {{1, 0, 0}, {0, Cos[O23], Sin[O23]}, {0, -Sin[O23], Cos[O23]}};
R2 = {{Cos[O13], 0, Sin[O13]*Exp[-I*dcp]}, {0, 1, 
    0}, {-Sin[O13]*Exp[I*dcp], 0, Cos[O13]}};
R3 = {{Cos[O12], Sin[O12], 0}, {-Sin[O12], Cos[O12] , 0}, {0, 0, 1}};
U = R1*R2*R3;
Ve[x_] = Sqrt[2]*Gf*Exp[-2*x/Xmax];
P[x_] = 1/(2 p) * U*M*
    ConjugateTranspose[U] + {{Ve[x], 0, 0}, {0, 0, 0}, {0, 0, 0}};
A[x_] = Table[ai[i, x], {i, 1, 3}];
A0 = Table[0, {3}];
A0[[1]] = 1;

diffEq = Table[
  A'[x][[i]] == -I (P[x].A[x])[[i]], {i, 1, 3}]; initialCond = 
 Table[A[0][[i]] == A0[[i]], {i, 1, 3}];
sol = NDSolve[{Join[diffEq, initialCond]}, A[x], {x, 0, Xmax}];
Plot[Evaluate[
  Table[(ai[i, x] /. sol)*Conjugate[ai[i, x] /. sol], {i, 1, 3}]], {x,
   0, Xmax}]

The solution is supposed to be

Thank you very much for your help !