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"];
  • $\begingroup$ First of all, you can already solve for R(t) by hand. Second, the remaining equation probably should have C replaced by C_i. Is that what you mean? You will then also need the initial conditions for all C_i, I would guess. With that, you would have a coupled system of first-order equations in time for the C_i. $\endgroup$
    – Jens
    Commented May 19, 2012 at 4:50
  • $\begingroup$ I've updated the question. Wb is a function dependent on C, I'm not sure if it possible to solve it by hand. $\endgroup$
    – Andrew
    Commented May 19, 2012 at 14:20
  • 1
    $\begingroup$ So, initially $x=x(y,z)$ but you would like to replace this 2d function by a finite number of 1d functions, ie, you'd like to discretise the $z$ direction by hand. Right? $\endgroup$
    – acl
    Commented May 19, 2012 at 17:07
  • 1
    $\begingroup$ @Andrew then you need to solve a set of differential equations for $x_i(y,t)$ (ie, a set of variables, not a single var) and $y(x_1,\ldots,x_N,t)$, whereas you're trying to set the problem up for a single $x$. $\endgroup$
    – acl
    Commented May 19, 2012 at 17:30
  • 1
    $\begingroup$ Look up Delayed Differential Equations in the documentation (howto/SolveDelayDifferentialEquations). Your ?? will be x[t-h]. Also figure out the correct initial condition because you will need to specify it for a range of values. $\endgroup$ Commented May 20, 2012 at 16:10

2 Answers 2


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.


(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:

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
        {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:

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

Mathematica graphics

Not sure if I am answering the right question though.

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

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

Mathematica graphics

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.

  • 1
    $\begingroup$ You know, you could have used EvaluationMonitor :> Sow[{t, y[t]}] as an NDSolve[] option setting, don't you? $\endgroup$ Commented May 19, 2012 at 18:32
  • 1
    $\begingroup$ @J.M. yes, but I am not the most elegant of programmers :) $\endgroup$
    – acl
    Commented May 19, 2012 at 18:35

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

enter image description here

  • $\begingroup$ Well, no, the problem is that I don't know what x should be. $\endgroup$
    – Andrew
    Commented May 19, 2012 at 14:18

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.