Clear["Global`*"];
    n = 30; ke = 0; W = 1.22; wh = 0.12; Wd = 0.005; xx = 4.125; uxing =
    0.028; h = 0.072;
    xm = xx/W;
    fai = Table[Cos[m*Pi*y], {m, 0, n}]
    NN = N[Integrate[fai^2, {y, 0, 1}]]
    a1 = 0.479;
    b1 = 0.03285;
    c1 = -0.1011;
    d1 = -31.6;
    uu1 = a1*Exp[b1*yy] + c1*Exp[d1*yy]
    a2 = -2.461*10^-16;
    b2 = 27.55;
    c2 = 0.4954;
    d2 = -0.02472;
    uu2 = a2*Exp[b2*yy] + c2*Exp[d2*yy]
    uav = Integrate[uu1, {yy, 0, W/2}]/W + Integrate[uu2, {yy, W/2, W}]/W
    e = 4*10^(-4);
    u1 = a1*Exp[b1*y*W] + c1*Exp[d1*y*W]
    u2 = a2*Exp[b2*y*W] + c2*Exp[d2*y*W]
    U1 = u1/uav
    U2 = u2/uav
    V = 0;
    dy = e/(uav*W); t = ke*W/uav; h1 = wh/W; h2 = (Wd + wh)/W;
    obsdatao =
      Import["https://pastebin.com/raw/Fz38738r", "Table"];
    obsdata = First[obsdatao];
    c0a = obsdata[[All, 1]]
    c0 = Sum[c0a[[i]]*0.01, {i, 1, 122}]/W
    eqns = Table[{FullSimplify[
         Sum[cc[j]'[x]/Sqrt[NN[[j + 1]]*NN[[m + 1]]]*
           Chop[
            Integrate[fai[[j + 1]]*fai[[m + 1]]*U1, {y, 0, W/2}] +
             Integrate[fai[[j + 1]]*fai[[m + 1]]*U2, {y, W/2, W}]], {j, 0,
            n}]] ==
        FullSimplify[
         Sum[cc[j][x]/
             Sqrt[NN[[j + 1]]*NN[[m + 1]]]*(Chop[
               Integrate[
                fai[[j + 1]]*D[dy*Dt[fai[[m + 1]], y], y], {y, 0, 1}]] -
              Integrate[
               V*fai[[m + 1]]*Dt[fai[[j + 1]], y], {y, 0, 1}]), {j, 0,
            n}] - t*cc[m][x]],
       cc[m][0] ==
        Sum[Integrate[
          fai[[m + 1]]/Sqrt[NN[[m + 1]]]*c0a[[k]]/c0, {y, (k - 1)/122,
           k/122}], {k, 1, 122}]}, {m, 0, n}]
    vars = Table[cc[j], {j, 0, n}];
    s = NDSolve[eqns, vars, {x, 0, xm}, PrecisionGoal -> Infinity,
       Method -> {(*"StiffnessSwitching",Automatic,*)
         "IndexReduction" -> Automatic,
         "EquationSimplification" -> {(*Automatic,*)"Residual",
           "TimeConstraint" -> 10, "SimplifySystem" -> True}}];
    s1 = Flatten[s]
   (* Table[cc[0][x] /. s1, {x, 0, xm}]
    Plot[cc[0][x] /. s1, {x, 0, xm}]
    Table[Plot[cc[m][x] /. s1, {x, 0, xm}], {m, 0, n}]
    Table[cc[m][400/122 + 1/2*25/122] /. s1, {m, 0, n}]*)
    c3 = Sum[fai[[m + 1]]/Sqrt[NN[[m + 1]]]*(Re[cc[m][400/122 + 1/2*25/122] /. s1]), {m, 0,
        n}]
    N[Table[c3, {y, 1/2*1/122, (1 - 1/2*1/122), 1/122}]]
    Plot[c3, {y, 1/2*1/122, (1 - 1/2*1/122)}, AxesLabel -> {y, c}]