I am currently trying to solve a time dependent Schrodinger equation (TDSE) for a rectangular potential barrier and a varying magnetic potential with a time dependent magnetic potential is oriented on the x-Axis and has a sinusoidal form.

The equation I'm trying to solve is: $$ i\omega\partial_{\xi}u(x,\xi) = -\frac{\hbar^2 c^2}{2\mu c^2}\partial_{x,x}u(x, \xi) + V(x)u(x,\xi) + i \frac{\hbar c}{\mu c^2}q* A0*cos(\xi)\partial_{x}u(x,\xi) $$

$\xi = \omega t$.

The initial state of the particle is a stationary state with energy Q > 0 found by:

hbarc = 197.326;
mu = 4000.;
S0 = 3.0;
omega = 0.0001;
ec = 0.0854;
q = 2.*ec;
Q = 5.0;
V0 = 0.;
V1 = 20.;
a = 10.;
l = 5.;
x0 = 0.;
xN = a + l + 10.;
A0 = mu omega S0/(hbarc q);
k = Sqrt[2 mu Q]/hbarc;
potential[x_] = Piecewise[{{V1,a<=x<=a+l},{V0,x<a||x>a+l}};

sol=ParametricNDSolve[Join[eqSchUnpert,bc],u0,{x, x0, xN},beta]

and, since it is not a bound state, any beta will do so I choose 10^-4:


Now, when adding the magnetic potential the code takes forever and doesn't bring any result:

eqSchPert = {I*omega*
 D[u[t, x],t]==-hbarc*hbarc/(2 mu)D[u[t, x],{x, 2}]+
 potential[x] u[t, x] + I*S0*omega*Cos[t]*D[u[t, x], x]};

with the initial and boundary contitions:

icPert = {u[0.,x]== psi0[x]};
bcPert = {u[t,x0+S0]== psi0[x0+S0 + S0*Sin[t]]*Exp[-I*Q*t/omega], u[t,xN-S0]== psi0[xN-S0 + S0*Sin[t]]*Exp[-I*Q*t/omega]};

NDSolve never finishes. If anyone knows what's wrong, please tell me. What I'm doing now is:

eqSyst=Join[eqSchPert, icPert, bcPert];

NDSolve[eqSyst,u,{t,0.,2.},{x,x0+S0,xN-S0},Method->{"MethodOfLines","SpatialDiscretization"->{"TensorProductGrid","DifferenceOrder"->"Pseudospectral","MinPoints"->300,"MaxPoints"-> 300}}]


After searching a little more I gave FiniteElement a try and I think I'm getting close, but I'm not there yet. If I set $\omega$ = $10^{4}$ and solve the equation with:

omega = 10.^4;
icPert = {u[r, 0.] == psi0[r]};
bcPert = {u[k*(x0 + S0), t] == 
psi0[k (x0 + S0 + S0*Sin[t])]*Exp[-I*Q*t/omega],
    u[k*(xN - S0), t] == 
psi0[k (xN - S0 + S0*Sin[t])]*Exp[-I*Q*t/omega]};
schEq = {I*omega*D[u[r, t], t] == -Q (D[u[r, t], {r, 2}]) + I*k*S0*omega*Cos[t] D[u[r, t], r] + potential[r]*u[r, t]};
eqSyst = Join[schEq, icPert, bcPert];

and I use a mesh made like this:

mesh = ToElementMesh[Rectangle[{k (x0 + S0), 0.}, {k (xN - S0), tFinal}], "MaxBoundaryCellMeasure" -> 0.05, "MeshElementType" -> QuadElement];

the solution becomes:

solution = NDSolve[eqSyst, u, {r, t} \[Element] mesh];
u1 = u /. First[solution];
Plot3D[Abs[Evaluate[u1[r, t]]]^2, {r, k (x0 + S0), k (xN - S0)}, {t, 0., 2*\[Pi]}, PlotRange -> All, Mesh -> All]

I obtain: enter image description here

which is exactly what I expected for $\omega$ very small. If I try this with the previous $\omega$ I get this:

enter image description here

which is absurd. In both cases, I receive a warning message:

NDSolve::femcscd: The PDE is convection dominated and the result may not be stable. Adding artificial diffusion may help.

Please note that in the equations I passed to the dimensionless variable $x->\rho = k x$, with $k$ defined in the beggining of this post.

  • $\begingroup$ Hi Stefan and welcome! To make the most of Mma.SE start by taking the tour now. It will help us to help you if you write an excellent question. Edit if improvable, show due diligence by giving a minimal working example of code you have tried yourself. As you receive give back, vote and answer questions, keep the site useful, be kind, correct mistakes and share what you have learned. $\endgroup$
    – rhermans
    Jul 16, 2018 at 12:55
  • 1
    $\begingroup$ Definition of eqSyst ??? $\endgroup$ Jul 16, 2018 at 13:11
  • $\begingroup$ This is mistake u[t,x0+S0]== psi0[x0+S0 + S0*Sin[t]]*Exp[-I*Q*t/omega]. $\endgroup$ Jul 16, 2018 at 13:44
  • $\begingroup$ Sorry, Ulrich, I just added the definition of eqSyst. $\endgroup$
    – Kvothe
    Jul 16, 2018 at 14:47
  • $\begingroup$ Alex Trounev, what do you mean? What's wrong? $\endgroup$
    – Kvothe
    Jul 16, 2018 at 14:49

1 Answer 1


Minimum working code, in which all equations have solutions. I made a replacement t/omega->t to reduce computation time.

hbarc = 197.326;
mu = 4000.;
S0 = 3.0;
omega = 0.0001;
ec = 0.0854;
q = 2.*ec;
Q = 5.0;
V0 = 0.;
V1 = 20.;
a = 10.;
l = 5.;
x0 = 0.;
xN = a + l + 10.;
A0 = mu omega S0/(hbarc q);
k = Sqrt[2 mu Q]/hbarc;
potential[x_] = 
  Piecewise[{{V1, a <= x <= a + l}, {V0, x < a || x > a + l}}];
eqSchUnpert = {Q*u0[x] == -hbarc*hbarc/(2*mu)*D[u0[x], {x, 2}] + 
bc = {u0[xN] == beta*Exp[I*k*xN], (D[u0[x], x] /. x -> xN) == 
psi0 = ParametricNDSolveValue[{eqSchUnpert, bc}, 
   u0, {x, x0, xN}, {beta}];
Plot[Abs[psi0[10^-4][x]], {x, x0, xN}]
eqSchPert = {I*
    D[u[t, x], t] == -hbarc*hbarc/(2 mu) D[u[t, x], {x, 2}] + 
    potential[x] u[t, x] + 
    I*S0*omega*Cos[t*omega]*D[u[t, x], x]}; icPert = {u[0., x] == 
bcPert = {u[t, x0 + S0] == psi0[10^-4][x0 + S0], 
   u[t, xN - S0] == psi0[10^-4][xN - S0]};

sol1 = NDSolveValue[{eqSchPert, icPert, bcPert}, 
   u[t, x], {t, 0., 2}, {x, x0 + S0, xN - S0}, 
   Method -> {"MethodOfLines", 
     "SpatialDiscretization" -> {"TensorProductGrid", 
       "MinPoints" -> 300, "MaxPoints" -> 1000, 
       "DifferenceOrder" -> "Pseudospectral"}}, MaxSteps -> 10^6];
Plot3D[Abs[sol1], {t, 0, 2}, {x, x0 + S0, xN - S0}, Mesh -> None, 
 ColorFunction -> Hue, AxesLabel -> {"t", "x"}]


  • $\begingroup$ Thank you very much! I know that NDSolve is able to find a solution for short time intervals. I want one the range $\omega$ t from 0 to 2 Pi. This is the reason I post this question. I don't know if this is possible. $\endgroup$
    – Kvothe
    Jul 17, 2018 at 6:33
  • $\begingroup$ I have modified the boundary condition for the second problem and now it does not give a warning for inconsistent boundary and initial conditions. $\endgroup$
    – Kvothe
    Jul 17, 2018 at 6:59
  • $\begingroup$ Thank you, Stefan. I made a calculation with these boundary conditions. To count for long times, we need to reduce the equation to a dimensionless form, determine the region $ 0\le x \le 2$, and lower the omega to $2\pi / 100$. In this case, there is a chance to count up to $t=100$. I launched such a task, but it has not yet been counted to the end, but it requires fast CPU and a memory of at least 10 GB. $\endgroup$ Jul 17, 2018 at 13:11
  • $\begingroup$ Of course we can scale the problem. But it must be scaled all together: the energy of the particle needs to be 4 orders of magnitude higher than $\omega$ and so on. You can choose any system of units as long as the quantities follow the ratio in my original post. $\endgroup$
    – Kvothe
    Jul 17, 2018 at 13:23
  • $\begingroup$ Why do not you do it yourself? Now the task looks like a freak. It is necessary to give it a form convenient for calculation. So I could not calculate the task for t=100. The computer did not have enough memory. Although it seems that the task is simple. But the Mathematica is hanging on it. $\endgroup$ Jul 17, 2018 at 15:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.