Needs["NDSolve`FEM`"]
    Domain = ImplicitRegion[
       0 <= x <= 1.25*10^-2 && 0 <= y <= 1*10^-2, {x, y}];
    meshA = ToElementMesh[Domain, MaxCellMeasure -> {"Length" -> 0.0007},
      "MaxBoundaryCellMeasure" -> 0.0002]
    meshA["Wireframe"]

    PDEtoMatrix[{pde_, Ga_}, u_, r_] :=
 Module[{ndstate, feData, sd, bcData, methodData,
   pdeData}, {ndstate} =
   NDSolve`ProcessEquations[Flatten[{pde, Ga}], u, Sequence @@ {r}];
  sd = ndstate["SolutionData"][[1]];
  feData = ndstate["FiniteElementData"];
  pdeData = feData["PDECoefficientData"];
  bcData = feData["BoundaryConditionData"];
  methodData = feData["FEMMethodData"];
  {DiscretizePDE[pdeData, methodData, sd],
   DiscretizeBoundaryConditions[bcData, methodData, sd], sd,
   methodData}]

Tin = 550;
pre = 10^4;
Ea = 5000*1000;
R = 1.986*1000;
k = pre*Exp[-Ea/Tin/R];
a = 0.5;
Ga = 3*10^-8;
fm = {.15, .10};
Rho = 1000;
PM0 = {24, 28};
CaZ = Rho*fm[[1]]/PM0[[1]];
CbZ = Rho*fm[[2]]/PM0[[2]];
Vx[x_, y_] := -a*(x - (1.25*10^-2)/2);
Vy[x_, y_] := a*(y - 0.005);
pde = {D[Ca[x, y]*Vx[x, y], x] + D[Ca[x, y]*Vy[x, y], y] -
     Ga*D[Ca [x, y], {x, 2}] - Ga*D[Ca [x, y], {y, 2}] == 0,
   D[Cb[x, y]*Vx[x, y], x] + D[Cb[x, y]*Vy[x, y], y] -
     Ga*D[Cb [x, y], {x, 2}] - Ga*D[Cb [x, y], {y, 2}] == 0};
bds = {DirichletCondition[Ca[x, y] == CaZ, x == 0],
   DirichletCondition[Ca[x, y] == 0, x == 1.25*10^-2],
   DirichletCondition[Cb[x, y] == CbZ, x == 1.25*10^-2],
   DirichletCondition[Cb[x, y] == 0, x == 0]};

{dPDE, dBC, sd, md} =
  PDEtoMatrix[{pde, bds}, {Ca, Cb}, {x, y} [Element] meshA,
   Method -> {"FiniteElement",
     "InterpolationOrder" -> {Ca -> 2, Cb -> 2},
     "MeshOptions" -> {"ImproveBoundaryPosition" -> False,
       "MaxCellMeasure" -> 0.001}}];

linearLoad = dPDE["LoadVector"];
linearStiffness = dPDE["StiffnessMatrix"];
vd = md["VariableData"];
offsets = md["IncidentOffsets"];

uOld = ConstantArray[{0.}, md["DegreesOfFreedom"]];
mesh2 = md["ElementMesh"];
mesh1 = MeshOrderAlteration[mesh2, 1];

ClearAll[rhs]
rhs[t_?NumericQ, ut_] := Module[{uOld}, uOld = ut;
  Do[ClearAll[Ca0, Cb0];

   (*create functions interpolations*)
   Ca0 = ElementMeshInterpolation[{mesh2},
     uOld[[offsets[[1]] + 1 ;; offsets[[2]]]]];
   Cb0 = ElementMeshInterpolation[{mesh2},
     uOld[[offsets[[2]] + 1 ;; offsets[[3]]]]];

   (*these are the linearized coefficients*)
   nlPdeCoeff =
    InitializePDECoefficients[vd, sd,
     "ConvectionCoefficients" -> {{{Vx[x, y], Vy[x, y]},{0,0}},{{0, 0}, {Vx[x, y], Vy[x, y]}}},
     "DiffusionCoefficients" -> {{{{-Ga, 0}, {0, -Ga}}, {{0, 0}, {0, 0}}}, {{{0, 0}, {0, 0}}, {{-Ga,
    0}, {0, -Ga}}}},
     "ReactionCoefficients" -> {{k*Cb0[x,y], 0}, {0, 0}}, {{0, 0},{0,k*Ca0[x,y]}}];

   nlsys = DiscretizePDE[nlPdeCoeff, md, sd];
   nlLoad = nlsys["LoadVector"];
   nlStiffness = nlsys["StiffnessMatrix"];
   ns = nlStiffness + linearStiffness;
   nl = nlLoad + linearLoad;
   DeployBoundaryConditions[{nl, ns}, dBC];
   diriPos = dBC["DirichletRows"];
   nl[[diriPos]] = nl[[diriPos]] - uOld[[diriPos]];

   dU = LinearSolve[N[ns], N[nl]];
   Print[i, " Residual: ", Norm[nl, Infinity], "  Correction: ",
    Norm[dU, Infinity]];
   uOld = uOld + dU;
   If[Norm[dU, Infinity] < 10^-6, Break[]];, {i, 8}];
      uOld]

uNew = rhs[0, uOld];