I would like to solve a 1-dimensional mass transport problem of the form:
$$\frac{\partial{c}}{\partial{t}}=-\nabla\cdot\mathbf{N}$$
where there are no homogeneous reactions and with initial and boundary conditions:
$$ c(t=0) = c_{bulk}$$
$$ c(y=0) = 0 $$
$$ c(y=\infty) = c_{bulk}$$
$c$ is concentration, $t$ is time, $\mathbf{N}$ is the flux, $y$ is the spatial coordinate, and $c_{bulk}$ is the bulk concentration
This problem can be re-written as an ordinary differential equation (ODE) using a similarity transform:
$$\frac{d\mathscr{N}}{d\eta}=2\eta\frac{dc}{d\eta}$$
where $\eta$ is defined as
$$\eta=\frac{y}{2\sqrt{Dt}}$$
where $D$ is the diffusion coefficient.
The relationship between $\mathscr{N}$ and $\mathbf{N}$ is given by
$$\mathscr{N}=2\sqrt{\frac{t}{D}}\mathbf{N}=-\frac{dc}{d\eta}$$
when the flux $\mathbf{N}$ follows Fickian diffusion
$$\mathbf{N}=-D\nabla{c}$$
The transformed boundary conditions are
$$c(\eta=0) = 0$$
$$c(\eta=\infty) = c_{bulk}$$
($\eta = 6$ is used to approximate infinity; note the numerical solution compares favorably with the analytical solution for $\eta=6$.)
substituting $\mathscr{N}=-dc/d\eta$ into the transformed mass balance equation produces
$$-\frac{d^2c}{d\eta^2}=2\eta\frac{dc}{d\eta}$$
which conveniently has an analytical solution for these boundary conditions: $c(\eta)=c_{\infty} \cdot \mathrm{erf} (\eta)$.
My question deals with solving this problem numerically as a boundary-value problem. Using central-differences to calculate the first and second spatial derivatives, i.e.
$$\frac{dc}{d\eta}=\frac{c_{j+1}-c_{j-1}}{2h}$$
$$\frac{d^2c}{d\eta^2}=\frac{c_{j+1}-2c_{j} + c_{j-1}}{h^2} $$
$ j = 1$ $$\frac{dc}{d\eta} = \frac{-3 c_j + 4 c_{j+1} - c_{j+2}}{2h}$$
$ j = NJ$ $$ \frac{dc}{d\eta} = \frac{3c_j - 4c_{j-1} + c_{j-2}}{2h} $$
I have no problem solving this equation: $-\frac{d^2c}{d\eta^2}=2\eta\frac{dc}{d\eta}$
However, I wish to segregate the equations and solve as a set of coupled first order ODEs:
$$\frac{d\mathscr{N}}{d\eta}=2\eta\frac{dc}{d\eta}$$
$$\mathscr{N}=-\frac{dc}{d\eta} $$
(2 first-order ODEs - still only need two boundary conditions), I get some oscillatory behavior in the solution which makes convergence difficult and sometimes practically impossible.
My question is what is causing this oscillatory behavior and how can I mitigate it?