# solving PDE with state-dependent boundary conditions

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 and change the boundary conditions, in odefun or odestep? I am asking this because the odefun routine might be used internally by the solver for evaluation of jacobian etc and ideally, we should not change boundary conditions in the middle of an integration step, right?

PS. I was not sure which was the better place to ask this question, programming stack exchange or maths stack exchange. I asked it here because its essentially a math problem.

• try instead the scientific computing SX: scicomp.stackexchange.com – user225318 Apr 8 '15 at 18:10
• @user225318 OK. Will do that. Thanks. – Sandeep Apr 8 '15 at 18:25