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 !