#Explanation for the Warning

Notice the description for the warning `NDSolve::mconly` is

>NDSolve::mconly: **For the method IDA**, only machine real code is available. 

What's method IDA? It's a method for solving DAE. (Please search in the document for more information. ) 

Why is `NDSolve` calling a DAE solver? Because `NDSolve` discretizes the PDE system based on method of lines (`"MethodOfLines"`) and the resulting system is a DAE system or an ODE system.

Why does `NDSolve` choose to discretize the PDE system to a DAE system rather than an ODE system? Because your PDE system involves terms like $\frac{\partial ^2v}{\partial t\partial y}$ and current implementation of method of lines in `NDSolve` isn't clever enough to transform the discretized system to the standard form required by the ODE solver. (For more information about the standard form, check [this post](https://mathematica.stackexchange.com/a/158519/1871). )

Why does the DAE solver fail? Because generally the DAE solver of _Mathematica_ is weaker than the ODE solver, at least up to now. 

#Solutions

##Partly `NDSolve`-based solution

One possible method to resolve the problem is to discretize the system ourselves and make `NDSolve` use the ODE solver as I've done [here](https://mathematica.stackexchange.com/a/163956/1871). I'll use [`pdetoode`](http://mathematica.stackexchange.com/a/127997/1871) for the task.

Notice I've modified the definition of `η0` and `v0` a little, because `pdetoode` can only handle `Listable` function.

    Clear[η0, v0, α, β]
    SetAttributes[#, Listable] & /@ {η0, v0};
    SeedRandom[1];
    η0[y_?NumericQ] = 
      BSplineFunction[Join[{0.}, RandomReal[{-1, 1}, 10], {0.}], 
        SplineClosed -> False][(y + 1)/2];
    SeedRandom[2];
    v0[y_?NumericQ] = 
      BSplineFunction[Join[{0.}, RandomReal[{-1, 1}, 10], {0.}], 
        SplineClosed -> False][(y + 1)/2];

    α = 1; β = 0.5; m = 300; Tend = 20; nx = 201;
    
    With[{η = η[t, y], v = v[t, y], U = 1 - y^2}, 
       feq = {D[η, t] + I α U η + op1[y, α, β][η]/m == (-I) β D[U, y] v, 
              op1[y, α, β][D[v, t]] + 
                I α U op1[y, α, β][v] + I α D[U, {y, 2}] v + op2[y, α, β][v]/m == 0}; 
     fic = {η == η0[y], v == v0[y]} /. t -> 0; 
     fbc = {{v == 0, η == 0, D[v, y] == 0} /. y -> -1, 
            {v == 0, η == 0, D[v, y] == 0} /. y -> 1}; ]

    domain = {-1, 1};
    difforder = 2;
    points = 50;
    grid = Array[# &, points, domain];
    (* Definition of pdetoode isn't included in this post,
       please find it in the link above. *)
    ptoofunc = pdetoode[{η, v}[t, y], t, grid, difforder];
    
    delone = #[[2 ;; -2]] &;
    deltwo = #[[3 ;; -3]] &;
    
    ode@1 = delone@ptoofunc@feq[[1]];
    ode@2 = deltwo@ptoofunc@feq[[2]];
    odeic = ptoofunc@fic;
    odebc = ptoofunc@With[{sf = 1}, Map[sf # + D[#, t] &, fbc, {3}]];
    var = Outer[#[#2] &, {η, v}, grid];
    
    sollst = NDSolveValue[{ode /@ {1, 2}, odeic, odebc}, var, {t, 0, Tend}];

    (* A more advanced but faster approach: *)
    (*
    lhs = D[Flatten[var][t] // Through, t];
    {barray, marray} = CoefficientArrays[{ode /@ {1, 2}, odebc} // Flatten, lhs];
    rhs = LinearSolve[marray, -barray]; // AbsoluteTiming
    
    sollst = NDSolveValue[{lhs == rhs, odeic}, var, {t, 0, Tend}];
     *)

    rebuild = ListInterpolation[
        Developer`ToPackedArray@#["ValuesOnGrid"] & /@ # // 
         Transpose, {#[[1]]["Coordinates"][[1]], grid}] &;
    
    fsol = rebuild /@ sollst;
    
    u1[t, y] = I/α D[v[t, y], y] /. v -> fsol[[2]];
    u3[t, y] = I/α η[t, y] /. η -> fsol[[1]];
    
    Plot3D[u1[t, y] // Abs // Evaluate, {t, 0, Tend}, {y, -1, 1}]

![Mathematica graphics](https://i.sstatic.net/ttgtv.png)

    Plot3D[u3[t, y] // Abs // Evaluate, {t, 0, Tend}, {y, -1, 1}]

![Mathematica graphics](https://i.sstatic.net/qnmUz.png)

You may try larger `difforder` or `points`, but notice this method will probably fail if their values are too large, because it'll be harder to transform the system to the standard form as they become larger.

If you want to make the solution fit the b.c. better, choose a larger `sf`. As to the meaning of `sf`, check [this post](https://mathematica.stackexchange.com/a/127411/1871).

##Purely FDM-based solution

As already mentioned, transforming the system to the standard form required by ODE solver can be time consuming. (`difforder = 2; points = 100` is already challenging for the method above. ) So it's not a bad ieda to leave `NDSolve` alone and turn to pure finite difference method (FDM) as I've done [here](https://mathematica.stackexchange.com/a/183767/1871). I'll use `pdetoae` for the task:

    (* Definitions for feq etc. are the same as above. *)
    Clear[domain, points, grid];
    domain@t = {0, Tend}; domain@y = {-1, 1};
    points@t = 100; points@y = 100;
    difforder = 4;
    (grid@# = Array[# &, points@#, domain@#]) & /@ {t, y};
    (* Definition of pdetoode isn't included in this post,
       please find it in the link above. *)
    ptoafunc = pdetoae[{η, v}[t, y], grid /@ {t, y}, difforder];
    delone = #[[2 ;; -2]] &;
    deltwo = #[[3 ;; -3]] &;
    ae@1 = delone /@ Rest@ptoafunc@feq[[1]];
    ae@2 = deltwo /@ Rest@ptoafunc@feq[[2]];
    aeic@1 = delone@ptoafunc@fic[[1]];
    aeic@2 = deltwo@ptoafunc@fic[[2]];
    aebc = ptoafunc@fbc;
    
    var = Outer[#[#2, #3] &, {η, v}, grid@t, grid@y];
    {barray, marray} = 
      CoefficientArrays[{Outer[#@#2 &, {ae, aeic}, {1, 2}], aebc} // Flatten, var // Flatten];
    sollst = LinearSolve[marray, -barray];
    solfunclst = 
      ListInterpolation[#, grid /@ {t, y}] & /@ ArrayReshape[sollst, {2, points@t, points@y}];
    u1[t, y] = I/α D[v[t, y], y] /. v -> solfunclst[[2]];
    u3[t, y] = I/α η[t, y] /. η -> solfunclst[[1]];
    
The resulting pictures look the same as above so I'd like to omit them here.