I am interested in solving the following PDE (heat equation): $$\frac{\partial u}{\partial t} = \kappa \frac{\partial ^2 u}{\partial x^2}$$ In order to solve it, I discretize space uniformly into $N$ segments and convert the PDE into $N+1$ ODEs: $$ \frac{du_k}{dt} = \kappa\frac{u_{k+1} + u_{k-1} - 2u_k}{(\Delta x)^2} \quad k=0,1\ldots.N$$ which I can solve using any ODE routine (such as ode45 or ode15s in MATLAB). For the particular case of my problem, the Neumann boundary conditions for the PDE are dependent on the values of $u_k$ itself, i.e $$ \frac{du_b}{dx}= \left\{ \begin{array}{cc} \alpha & u_b < u_{critical}, \\ -\gamma u_b & u_b \geq u_{critical} \end{array} \right. $$ where $\alpha$ and $\gamma$ are constants and $b=\{0,N\}$.

Assuming that I am solving the above problem using some numerical ODE solver (such as ode45 in MATLAB), and there are two routines: odefun (which I provide to the solver) to evaluate the derivatives and odestep which I define and the solver calls it after every successful integration step, where should I check if the boundary condition needs to be changed, in odefun or odestep? I am asking this because the odefun routine might be used internally by the solver for evaluation of jacobian, errors etc and ideally, we should not change boundary conditions in the middle of an integration step, right?

  • $\begingroup$ Welcome to StackExchange! I think this question is off-topic here because it is very specific to Matlab, not a general question about computational science. I think this question would be better placed on a Matlab specific forum. $\endgroup$ Commented Apr 10, 2015 at 12:38
  • $\begingroup$ @WolfgangBangerth That's partly true. However, my main question is about how to check for change in boundary conditions (for any numerical solver), and rather than using an abstract solver, I used Matlab's routines as an example. I was confused whether to ask it here or not, but just decided to give it a try. $\endgroup$
    – Chatter
    Commented Apr 10, 2015 at 16:37
  • $\begingroup$ I don't know anything about the odestep function, so I can't answer your question as posed. If you rephrased it as, say, how to implement the solution using Euler's method, or RK4, or some integrator you implement yourself, then I could probably help you. But you should also specify how you are going to impose the boundary conditions numerically. $\endgroup$ Commented Apr 11, 2015 at 3:55

1 Answer 1


You shouldn't have any problem coding in a conditional boundary condition into the odefun, like the one you defined in your problem statement.

Just use the last timestep's value of $u_b$ to do the conditional statement and grab the correct boundary condition value.

And the reality is, if your problem needs to change the boundary condition in the middle of integrating, then it has to do it. That's just how the physics are being modeled in your case. Given that these dynamics are okay, you should ideally incorporate that boundary condition change whenever you need to.


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.