(* Tested in v9.0.1 *)
Clear[η0, v0, α, β]

op1[y_, α_, β_] = ((α^2 + β^2)*# - D[#, {y, 2}]) &;
op2[y_, α_, β_] = (op1[y, α, β]@ op1[y, α, β]@#) &;

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}, diffbc[t, sf]@fbc];
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}];
 *)

fsol = rebuild[#, grid]& /@ sollst;

Clear[u1, u3]    
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}]

Remark

There seems to be a backslide of NDSolve at least since v12.3. I'm not sure what the root of the problem is, but if you find the code above doesn't produce the shown figure, use the following defintion of sollst

sollst = 
   NDSolveValue[Rationalize[#, 0] &@{ode /@ {1, 2}, odeic, odebc},
     var, {t, 0, Tend},
     SolveDelayed -> False, WorkingPrecision -> 16]; // AbsoluteTiming
(* {6.64703, Null} *)

or the more advanced but faster approach in the comment instead.

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.

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 idea to leave NDSolve alone and turn to pure finite difference method (FDM) as I've done here. I'll use pdetoae for the task:

(* Tested in v9.0.1 *)
Clear[η0, v0, α, β]

op1[y_, α_, β_] = ((α^2 + β^2)*# - D[#, {y, 2}]) &;
op2[y_, α_, β_] = (op1[y, α, β]@ op1[y, α, β]@#) &;

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}, diffbc[t, sf]@fbc];
var = Outer[#[#2] &, {η, v}, grid];

sollst = NDSolveValue[{ode /@ {1, 2}, odeic, odebc}, var, {t, 0, Tend}, 
  Method -> {"EquationSimplification" -> "MassMatrix"}];

fsol = rebuild[#, grid]& /@ sollst;

Clear[u1, u3]    
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}]

You may try larger difforder or points.

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 MassMatrix method shown above turns out to be very efficient for the problem. Still, it's not a bad idea to leave NDSolve alone and turn to pure finite difference method (FDM) as I've done here. I'll use pdetoae for the task:

Remark

There seems to be a backslide of NDSolve at least since v12.3. I'm not sure what the root of the problem is, but if you find the code above doesn't produce the shown figure, use the following defintion of sollst

sollst = 
   NDSolveValue[Rationalize[#, 0] &@{ode /@ {1, 2}, odeic, odebc},
     var, {t, 0, Tend},
     SolveDelayed -> False, WorkingPrecision -> 16]; // AbsoluteTiming
(* {6.64703, Null} *)

or the more advanced but faster approach in the comment instead.

#Explanation for the Warning

#Solutions

##Partly NDSolve-based solution

##Purely FDM-based solution

Explanation for the Warning

Solutions

Partly NDSolve-based solution

Purely FDM-based solution