# Using DSolve for a coupled differential equation

I'm trying to solve (using the DSolve command) the following system

$\psi_1''[x]+(E_1-E-xV)\psi_1[x]-i\psi_2'[x]=0$

$\psi_2''[x]+(E_2-E-xV)\psi_2[x]-i\psi_1'[x]=0$

However, seems the Mathematica cannot manage to give a solution. Am I doing something wrong? Do you guys have any suggestions concerning some kind of transformation which will give a simpler equation to be solved?

DSolve[{ψ1''[x] + (-E1 - En - x V0) ψ1[x] - I ψ2'[x] == 0,
ψ2''[x] + (-E2 - En - x V0) ψ2[x] - I ψ1'[x] == 0},  {ψ1, ψ2}, x]


Thanks too much!

Best,

• Please post the actual Mathematica code that you entered, properly formatted in code blocks. People like to copy and paste the code from the post into their own copies of Mathematica. Do you expect an analytic solution here (it looks possible: possibly Airy functions)? – march Dec 1 '16 at 22:08
• Schrödinger equations require boundary conditions to be fixed. So what are the conditions? – Alexei Boulbitch Dec 2 '16 at 9:21
• Thanks for the answer. Actually I still do not know how to input Mathematica codes properly.. I'll google it.. Yes, I'm expecting analytical solutions once I have the Airy functions considering $P=0$ Concerning the boundary conditions, I want $\psi1(0)=\psi2(0)=\psi1(L)=\psi1(L)=0$ – denis Dec 2 '16 at 12:06
• Are you confident that the ODE system has a symbolic solution? – bbgodfrey Dec 2 '16 at 17:22
• No so sure about that.. – denis Dec 5 '16 at 16:10

Actually it is not true that DSolve cannot solve it: just consider the uncoupled system formed by the sum $s[t]$ and the difference $d[t]$ of $\Psi_1[t]$ and $\Psi_2[t]$:

DSolve[{s''[t] + (e1 + e2 - 2 e - t v) s[t] - I s'[t] == 0,
d''[t] + (e1 - e2) d[t] + I d'[t] == 0}, {s, d}, t]


with solution:

(*{{s -> Function[{t},
E^((I t)/2)
AiryAi[(1/4 (-1 + 8 e - 4 e1 - 4 e2) + t v)/v^(2/3)] C +
E^((I t)/2)
AiryBi[(1/4 (-1 + 8 e - 4 e1 - 4 e2) + t v)/v^(2/3)] C],
d -> Function[{t},
E^(1/2 (-I - Sqrt[-1 - 4 e1 + 4 e2]) t) C +
E^(1/2 (-I + Sqrt[-1 - 4 e1 + 4 e2]) t) C]}}*)


Now you can get the solution for the original functions.

One way to find a solution easily is to go for numerical solution using NDSolve.

Here, I have choose random values for the parameters.

E1 = E2 = En = V0 = 1; L1 = 5;
sol = NDSolve[{ψ1''[x] + (-E1 - En - x V0) ψ1[x] - I ψ2'[x] == 0,
ψ2''[x] + (-E2 - En - x V0) ψ2[x] - I ψ1'[x] == 0, ψ1 == 0, ψ1[L1] == 0, ψ2 == 0,
ψ2[L1] == 0}, {ψ1, ψ2}, {x, 0, L1}];


Finally, plotting the results,

Plot[Evaluate[{Re[ψ1[x]], Re[ψ2[x]], Im[ψ1[x]],Im[ψ2[x]]} /. sol], {x, 0, L1},
PlotRange -> All, Frame -> True,PlotStyle -> {Red, Green, Directive[Blue, Dashed],
Directive[Black, Dashed]}] 