NDSolve and results of the previous computations

Question

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"];

Answer 1

It seems (from the comments) that what you want to do is this: initially x=x(y,z) but you would like to replace this function of 2 vars by a finite number of functions of 1 var, ie, you'd like to discretize the $z$ direction by hand. This means that you need to solve a set of differential equations for the $x_i(y,t)$ (a set of functions, not a single fun) and $y(x_1,\ldots,x_N,t)$.

You're trying to set the problem up for a single $x$, and that is the problem.

OLD ANSWER

(I think this is neat so I'll leave it here for now)

I very likely misunderstood you. If, however, I understood correctly the question, it is something like: Suppose I have $y'(t)=f(y)$ and ask mathematica to solve it numerically. How do I inspect the values of $t$ used?

The answer is to rig the ODE so that you can peek at the values mathematica evaluates, as follows:

ClearAll[f];
upT = 2*Pi;
f[y_?NumericQ, t_?NumericQ] := (Sow[{t, y}]; -Sin@y)

Here I define $f(y)=-\sin(y)$, but in a way that allows me to watch which values are passed to $f$ by NDSolve (the Sow bit). Then:

points = (sol = NDSolve[
        {
         y'[t] == f[y[t], t],
         y[0] \[Equal] 1
         },
        y,
        {t, 0, upT}
        ]) // Reap // Last // Last;

solves the ODE, collecting the values of $y$ and $t$ passed to $f$ by NDSolve. Then plot the solution, along with the points at which $f$ was evaluated:

Plot[
 y[t] /. sol,
 {t, 0, upT},
 Epilog :> {Red, PointSize[.015], Point[points]},
 ImageSize -> 640
 ]

Not sure if I am answering the right question though.

It's also interesting to look at the size of the steps taken:

ListPlot[
 Differences[points[[All, 1]]],
 ImageSize -> 640,
 BaseStyle -> FontSize -> 20,
 AxesLabel -> {"i", "t[i]-t[i-1]"}
 ]

so they are not all equal (not all $h$), sometimes the solver goes backwards etc. Of course this will depend on the solver you use.

Answer 2

Your syntax has me a bit confused. SetDelayed (:=) looks off to me. Typically one uses this to define functions. I think you want Equal (==) where you have the system of equations and Set (=) where you assign the value of the expression to "solution".

Does this get you any closer (I replaced your "??" with "x")?

solution = First@NDSolve[{
c'[t] == (2*Wb[c[t]])/(r[t]*Sin[theta])*(1 - C/gamma) + 
  v/(Pi*r[t]^2)*(c[t] - x)/h, 
r'[t] == Wb[c[t]]/(gamma*Sin[theta]),
c[0] == 0,
r[0] == 0.15
}, {c, r}, {t, 0, 25900000}, Method -> "BDF"]