I want to solve following system of ODEs:
$ \Bigg\{ \begin{array}{} \frac{\partial C}{\partial t}=\frac{2W_b}{Rsin(\theta)}(1-\frac{C}{\gamma})+\frac{\nu}{\pi R^2}\frac{\partial C}{\partial Z} \\ \frac{\partial R}{\partial t}=\frac{W_b}{\gamma sin(\theta)} \end{array} $
Where $W_b=k_bC_h(1-\frac{C}{C_h})^{\frac{4}{3}}$ and $\nu, k_b, C_h, h, \gamma, \theta$ are constants.
It is known that:
$ \frac{\partial C}{\partial Z}=\frac{C_i-C_{i-1}}{Z_i-Z_{i-1}} $, where $Z_i-Z_{i-1}=h$ and $h$ is a constant.
So, the first equation become this: $ \frac{\partial C}{\partial t}=\frac{2W_b}{Rsin(\theta)}(1-\frac{C}{\gamma})+\frac{\nu}{\pi R^2}\frac{C_i-C_{i-1}}{h} $
My question is how can I get the result of previous computation, namely $C_{i-1}$?
Here is the mathematica code I wrote:
Wb[C_] := kb*Ch*(1 - C/Ch)^(4/3);
system := {
c'[t] == (2*Wb[c[t]])/(r[t]*Sin[theta])*(1 - C/gamma) + v/(Pi*r[t]^2)*(c[t] - ?? )/h,
r'[t] == Wb[c[t]]/(gamma*Sin[theta]),
c[0] == 0, r[0] == 0.15};
solution = First@NDSolve[system, {c, r}, {t, 0, 25900000}, Method -> "BDF"];
UPDATE: Simplified version of the problem.
This is a model based on a plug flow reactor.
$
\Bigg\{
\begin{array}{}
\dot{x} = \frac{c_1 x f(x)}{y} + \frac{\partial x}{\partial z} \\
\dot{y} = c_2 f(x)
\end{array}
$
Where c1 and c2 are constants. This is a model of a physical process and it was shown that z
dimension is quantified by a chunks of a constant size h
and x
is monotonously increasing, so $\frac{\partial x}{\partial z}=\frac{\Delta x}{\Delta z}=\frac{x_i - x_{i-1}}{h}$
f[x_] := ...;
system := {
x'[t] == c1*x*f[x[t]]/y[t] + (x[t] - ?)/h,
y'[t] == c2*f[x[t]],
x[0] == 0, y[0] == 0.15};
solution = First@NDSolve[system, {x, y}, {t, 0, 25900000}, Method -> "BDF"];