PDEPE for simple coupled ODE-PDE system?

Hi all, I have a system of coupled reaction-diffusion ODEs and PDEs that I wish to use pdepe to solve, and have been running into unexpected outputs. To check, I used a simpler "test" system which could be solved in an exact fashion, purely to test the output of PDEPE.

It is simply:

 function [cc,ff,ss]=brpde(x,t,u,dstatedx)

 mu=0.02;

 u1=u(1);
 u2=u(2);

     if (x >5)
        du1 = -0.07*u1;
        else
            du1=0;
        end

        if (x > 5)
        du2=-0.07*u2;
        else
            du2=0;
        end

 ss=[du1; 
    du2];

 ff=[0;
    mu*dstatedx(2);];

 cc=[1;
    1;];

 return

I have set the boundary conditions to impose no flux and the left and right hand boundaries (this should work, I believe, even for the ODE variable u1 where flux is zero throughout the slab).

function [pa,qa,pb,qb]=brbc(xa,ua,xb,ub,t)

vecone=ones(2,1);
veczero=zeros(2,1);

   pa =veczero;
   qa = vecone; 
   pb = veczero;
   qb = vecone;
  return

And suitable initial conditions:

function u0=bric(x)
u0=[10;
   10;];
return

The output for u2 looks correct as I would expect it. However, although the exact solution is that u1 (the variable with no flux) should asymptotically approach zero for x>5 and stay at 10 for x<5, at the 'interface', ie near x=5, the solutions transiently dip below zero (even when the mesh is made finer or the error tolerances are made more stringent). It's confusing behavior, and I am wondering whether it's since pdepe is not suited for mixed ODE-PDE systems. When I code this simple system into a discretized forward Euler, of course the output is as I would predict it.

Thought I would reach out to the experience of the community. Anyone have suggestions? Or is there an obvious error in the above? (My real application has coupling between the ODEs and PDEs, but if it doesn't work for this test case, I am leaning toward hard-coding it).

Thanks in advance! VI