Strange wrong result for (un)coupled PDEs using MATLAB's pdepe, time is doubled

Question

I am trying to solve two coupled reaction diffusion equations in 1d, using pdpe, namely

$\partial_t u_1 = \nabla^2 u_1 + 2k(-u_1^2+u_2)$

$\partial_t u_2 = \nabla^2 u_1 + k(u_1^2-u_2)$

The solution is in the domain $x\in[0,1]$, with initial conditions being two identical Gaussian profiles centered at $x=1/2$. The boundary conditions are absorbing for both components, i.e. $u_1(0)=u_2(0)=u_1(1)=u_2(1)=0$.

Pdepe gives me a solution without prompting any errors. However, I think the solutions must be wrong, because when I set the coupling to zero, i.e. $k=0$ (and also if I set it to be very small, say $k=0.001$), the solutions do not coincide with the solution of the simple diffusion equation

$\partial_t u = \nabla^2 u$

as obtained from pdepe itself.

Strangely enough, the solutions $u_1(t)=u_2(t)$ from the "coupled" case with coupling set to zero, and the solution for the case uncoupled by construction $u(t')$ coincide if we set $t'=2t$, that is, the solution of the "coupled" case evolves twice as fast as the solution of the uncoupled case.

Here's a minimal working example:

Coupled case

    function [xmesh,tspan,sol] = coupled(k) %argument is the coupling k

    std=0.001; %width of initial gaussian
    center=1/2; %center of gaussian

    xmesh=linspace(0,1,10000);
    tspan=linspace(0,1,1000);

    sol = pdepe(0,@pdefun,@icfun,@bcfun,xmesh,tspan);

        function [c,f,s] = pdefun(x,t,u,dudx)
            c=ones(2,1);

            f=zeros(2,1);
            f(1) = dudx(1);
            f(2) = dudx(2);

            s=zeros(2,1);
            s(1) = 2*k*(u(2)-u(1)^2);
            s(2) = k*(u(1)^2-u(2));
        end

        function u0 = icfun(x)
            u0=ones(2,1);
            u0(1) = exp(-(x-center)^2/(2*std^2))/(sqrt(2*pi)*std);
            u0(2) = exp(-(x-center)^2/(2*std^2))/(sqrt(2*pi)*std);
        end

        function [pL,qL,pR,qR] = bcfun(xL,uL,xR,uR,t)
            pL=zeros(2,1);
            pL(1) = uL(1);
            pL(2) = uL(2);

            pR=zeros(2,1);
            pR(1) = uR(1);
            pR(2) = uR(2);

            qL = [0 0;0 0];
            qR = [0 0;0 0];
        end

    end

Uncoupled case

    function [xmesh,tspan,sol] = uncoupled()

    std=0.001; %width of initial gaussian
    center=1/2; %center of gaussian

    xmesh=linspace(0,1,10000);
    tspan=linspace(0,1,1000);

    sol = pdepe(0,@pdefun,@icfun,@bcfun,xmesh,tspan);

        function [c,f,s] = pdefun(x,t,u,dudx)
            c=1;

            f = dudx;

            s=0;
        end

        function u0 = icfun(x)
            u0=exp(-(x-center)^2/(2*std^2))/(sqrt(2*pi)*std);
        end

        function [pL,qL,pR,qR] = bcfun(xL,uL,xR,uR,t)
            pL=uL;

            pR=uR;

            qL = 0;
            qR = 0;
        end

    end

Now, suppose we run

    [xmesh,tspan,soluncoupled] = uncoupled();
    [xmesh,tspan,solcoupled] = coupled(0); %coupling k=0, i.e. uncoupled solutions

One can directly check by plotting the solutions for any time index $it$ that, even if they should be identical, the solutions given by each function are not identical, e.g.

    hold all
    plot(xmesh,soluncoupled(it+1,:),'b')
    plot(xmesh,solcoupled(it+1,:,1),'r')
    plot(xmesh,solcoupled(it+1,:,2),'g')

On the other hand, if we double the time of the uncoupled solution, the solutions are identical

    hold all
    plot(xmesh,soluncoupled(2*it+1,:),'b')
    plot(xmesh,solcoupled(it+1,:,1),'r')
    plot(xmesh,solcoupled(it+1,:,2),'g')

The case $k=0$ is not singular, one can set $k$ to be small but finite, and the deviations from the case $k=0$ are minimal, i.e. the solution still goes twice as fast as the uncoupled solution.

I really don't understand what is going on. I need to work on the coupled case, but obviously I don't trust the results if it does not give the right limit when $k\to 0$. I don't see where I could be making a mistake. Could it be a bug?

Answer 1

I found the source of the error. The problem lies in the qL and qR variables of bcfun for the coupled() function. The MATLAB documentation, see here and here, is slightly ambiguous on whether the q's should be matrices or column vectors. I had used matrices

        qL = [0 0;0 0];
        qR = [0 0;0 0];

but in reality I should have used column vectors

        qL = [0;0];
        qR = [0;0];

Amazingly, pdpe didn't throw an error, and simply gave wrong results. This should perhaps be fixed by the developers.