Why is ODE15s ignoring my boundary conditions for the Fishers equation?

Question

I'm attempting to solve the Fishers PDE on $(x,t) \in (0,1)\times (0,T)$ by spatially discretising (method of lines) and then feeding the time derivative to ODE15s (therefore creating a system of ODEs). Here's the equation in full;

$$\frac{\partial u}{\partial t} = \frac{\partial ^2 u}{\partial x^2} +u(1-u), \quad \quad \frac{\partial u}{\partial x}(0,t) = 0, \quad \quad u(1,t) = 0.$$ I've set the initial conditions as $n(x,0) = 0.6(1-x^2)$ and then I have set the second element of this to $0.6$ so it respects the boundary conditions at $x=0$. Here's my code:

space_steps=500;
T=100;
x = linspace(0,1,space_steps).';
T_span=[0;T]; %T_vector
n_0 = 0.6*(1-x.^2);
n_0(2)=0.6;
[T_out,Y_out] = ode15s(@(t,y) Fisher(t,y,space_steps),T_span,n_0);

Here's the function 'Fisher' that does all the work:

function output = Fisher(~,input,space_steps)
n = [input(2);input(2:end-1);0];    %dndx = 0 at x=0 so input(1)=input(2), n=0 at the x=1.
dx=1/(space_steps-1);               %Delta x
np = [n(2:end);NaN]; %n_{i+1}
nm=[NaN;n(1:end-1)]; %n_{i-1}
%Second order derivative with a forward and backward finite difference quotient at each boundary respectively
d2ndx2 = [(2*n(1)-5*n(2)+4*n(3)-n(4))/dx; np(2:end-1)-2*n(2:end-1)+nm(2:end-1);(2*n(end)-5*n(end-1)+4*n(end-2)-n(end-3))/dx]/dx^2;
dndt =  d2ndx2 + n.*(1-n); %sets up the equation 
output = dndt; %gives it to ODE15s
end 

Now for some reason the function or ODE15s is completely ignoring the two boundary conditions I am prescribing it at each run at the begining. Could anyone explain why this is and provide a suitable remedy? Thanks. I'd also appreciate any feedback on my code that you could offer.

Answer 1

You want to use $u_0(t)=u_1(t)$, $u_N(t)=0$ to implement the boundary conditions. That leaves the inner nodes as unknown state variables $u_1(t),...,u_{N-1}(t)$, so it will not be the solver that ignores the boundary values, you construct the code to ignore them from the outset in the computation and only reconstruct and insert them as intermediate values where they are needed. That is, in the computation of the time derivative for the inner nodes via the central second order difference quotient. There is then also no need for one-sided second-order difference quotients, as the differential equation is not evaluated for the outer nodes. You can do this with the code

x = linspace(0,1,space_steps+1).';
dx = x(2)-x(1);
g = @(x) 0.6*(1-x.^2);
n0 = g(x(2:end-1));

function dndt = Fisher(~, u, dx)
    f = n*(1-n);
    dndt = conv([ n(1); n; 0], [1,-2,1], shape="valid")/dx^2 + f;
end;