# Problem with NDSolve and NDSolve::ivone error

I was trying to solve an diffusion-reaction rate model for the oxidation of Silicon (based on the nice paper by Almeida et.al.) with the help of NDSolve. In this model the boundary conditions are a bit unusual compared to the way Mathematica expects them. The command I tried looks like the following:

res = NDSolve[{
D[nO2[x, t], {t, 1}] == D0 Exp[-(ED/(kB T))] D[nO2[x, t], {x, 2}] -
K0 Exp[-(EK/(kB T))] nO2[x, t] nSi[x, t],
nO2[0.0, t] == nO20[1, T],
D[nSi[x, t], {t, 1}] == -K0 E^(-(EK/(kB T))) nO2[x, t] nSi[x, t],
nSi[x, 0.0] == 1.0
} //. {
D0 -> 1.291 10^11, ED -> 2.22,
K0 -> 2.022 10^7, EK -> 1.42,
nSi0 -> 0.91 10^15,
kB -> QuantityMagnitude[UnitConvert[Quantity["BoltzmannConstant"/"ElementaryCharge"]]],
T -> 700 + 273.15
}, {nO2, nSi}, {x, 0, 3000}, {t, 0, 60 60}][[1]]


where Mathematica complains about the way how the boundary condition for the function nO2[x,t] has been defined and drops an NDSolve::ivone error.

Playing around with this boundary condition I came up with a setup which gave solutions to the functions nO2[x,t] and nSi[x,t], however the choice of the x coordinate for the upper boundary where nO2 is zero is causing a strong influence on the solution as one can see when NDSolve is run for two different settings for the x coordinate:

res = NDSolve[{
D[nO2[x, t], {t, 1}] == D0 Exp[-(ED/(kB T))] D[nO2[x, t], {x, 2}] -
K0 Exp[-(EK/(kB T))] nO2[x, t] nSi[x, t],
nO2[x, 0.0] == Piecewise[{{0.161/T, x == 0}}, 0],
nO2[3000.0, t] == 0.0,
nO2[0.0, t] == nO20[1, T],
D[nSi[x, t], {t, 1}] == -K0 E^(-(EK/(kB T))) nO2[x, t] nSi[x, t],
nSi[x, 0.0] == 1.0
} //. {
D0 -> 1.291 10^11, ED -> 2.22,
K0 -> 2.022 10^7, EK -> 1.42,
nSi0 -> 0.91 10^15,
kB -> QuantityMagnitude[UnitConvert[Quantity["BoltzmannConstant"/"ElementaryCharge"]]],
T -> 700 + 273.15
}, {nO2, nSi}, {x, 0, 3000}, {t, 0, 60 60}][[1]]
Plot[1 - (nSi /. res)[x, 60 60], {x, 0, 30}, PlotRange -> {0, 1.1}]


yields the following Graph:

whereas the following code

res = NDSolve[{
D[nO2[x, t], {t, 1}] == D0 Exp[-(ED/(kB T))] D[nO2[x, t], {x, 2}] -
K0 Exp[-(EK/(kB T))] nO2[x, t] nSi[x, t],
nO2[x, 0.0] == Piecewise[{{0.161/T, x == 0}}, 0],
nO2[300.0, t] == 0.0,
nO2[0.0, t] == nO20[1, T],
D[nSi[x, t], {t, 1}] == -K0 E^(-(EK/(kB T))) nO2[x, t] nSi[x, t],
nSi[x, 0.0] == 1.0
} //. {
D0 -> 1.291 10^11, ED -> 2.22,
K0 -> 2.022 10^7, EK -> 1.42,
nSi0 -> 0.91 10^15,
kB -> QuantityMagnitude[UnitConvert[Quantity["BoltzmannConstant"/"ElementaryCharge"]]],
T -> 700 + 273.15
}, {nO2, nSi}, {x, 0, 300}, {t, 0, 60 60}][[1]]
Plot[1 - (nSi /. res)[x, 60 60], {x, 0, 30}, PlotRange -> {0, 1.1}]


gives the following graph

It's obvious that one gets different results by moving the position of the boundary condition and thus using the boundary condition in this way is not an option. I would be interested in finding a possibility to specify the boundary condition for nO2 as shown in the first code snippet.

-