PDEPE: Spatial Discretization has failed

I want to solve this system of partial differential equations:

Subject to the following initial conditions:

And the following boundary conditions: BC 1 (at x=0): (Yes it is the same as the initial condition)

BC 2 (at x=delta):

And this is my code:

if true

 function pdesolveme(electrolyteconc,currentdens,delta)

 [AdjDiffCoeff,diffCoeff,initialEQVal,product,rateCoeff]=DataScript;

 %These just assign appropriate values to the following constants

 DCO2=AdjDiffCoeff{electrolyteconc,'CO2'};
 DHCO3=AdjDiffCoeff{electrolyteconc,'HCO3'};
 DCO3=AdjDiffCoeff{electrolyteconc,'CO3'};
 DOH=AdjDiffCoeff{electrolyteconc,'OH'};

 k1f=rateCoeff{'3b','fwd'};
 k1r=rateCoeff{'3b','rev'};
 k2f=rateCoeff{'4','fwd'};
 k2r=rateCoeff{'4','rev'};

 %PDE Solver

 m=0;
 x=linspace(0,delta,100);
 t=linspace(0,20,100');

 sol=pdepe(m,@myPDE,@myPDEic,@myPDEbc,x,t);
 CO2=sol(:,:,1);
 HCO3=sol(:,:,2);
 CO3=sol(:,:,3);
 OH=sol(:,:,4);

 %Plots

 pOH=-log10(OH);
 pH=14-pOH;

 figure
 surf(x,t,pH)
 title('pH')
 xlabel('Distance x')
 ylabel('Time t')

    function [ c,f,s ] = myPDE( x,t,u,dudx )

        c=[1;1;1;1];
        f=[DCO2*dudx(1);DHCO3*dudx(2);DCO3*dudx(3);DOH*dudx(4)];
        s=[-u(1)*u(4)*k1f+u(2)*k1r;...
            u(1)*u(4)*k1f-u(2)*k1r-u(2)*u(4)*k2f+u(3)*k2r;...
            u(2)*u(4)*k2f-u(3)*k2r;...
            -u(1)*u(4)*k1f+u(2)*k1r-u(2)*u(4)*k2f+u(3)*k2r];
    end

    function u0 = myPDEic (x)
        u0 = initialEQVal{electrolyteconc,1:4}';
    end

    function [p1,q1,pr,qr] = myPDEbc(xl,ul,xr,ur,t)
        CO2Cons=CO2Consumption(currentdens);
        OHForm=OHFormation(currentdens);

        p1=initialEQVal{electrolyteconc,1:4}';
        q1=[0;0;0;0];

        pr=[CO2Cons;0;0;-OHForm];
        qr=[1;1;1;1];
    end
end
  end

Upon running this code, I get the following error message:

Error using pdepe (line 293) Spatial discretization has failed. Discretization supports only parabolic and elliptic equations, with flux term involving spatial derivative.

Error in pdesolveme (line 22) sol=pdepe(m,@myPDE,@myPDEic,@myPDEbc,x,t);

Could you explain to me what I did wrong in my code?

[Torsten]

p1=initialEQVal{electrolyteconc,1:4}'-ul(1:4);

(Don't know whether " transpose " is correct or not in this case).

Best wishes

Torsten.