Nonlinear PDE with a dynamical parameter

I am trying to solve fourth order partial differential equation describing a fluid problem, by using Mathematica. The PDE is provided here as:

$V_{\eta\eta t}=\frac{1}{H}(2\dot H V_{\eta\eta} +\eta\dot H V_{\eta\eta\eta}) + \frac{1}{RH^2}V_{\eta\eta\eta\eta}+\frac{1}{H}(VV_{\eta\eta\eta} - V_\eta V_{\eta\eta})$

Along with the boundary conditions:

$V(1,t)=-\dot H,$ $V(-1,t)=\dot H,$ $V_\eta(1,t)=V_\eta (-1,t)=0$

Initial condition is not provided but they have choose the initial condition as mentioned in the article ${Page. No: 61}$ in the paragraph after equation $(2.6)$ and at the end of $Page. No: 61$.

$H=1+\Delta \cos(2t)$ and $\dot{H}$ is the derivative of $H$ with respect to $t$.

My working for the given problem is:

    Attributes[MakeVariables] = {Listable};
MakeVariables[var_, n_] := Table[Unique[var], {n}]

soln[delta_] := 
 Module[{\[CapitalDelta] = delta, dt = 1/50, nmax = 50, dy = 1/10}, 
  H[t_] := 1 + \[CapitalDelta]*Cos[2*t];
  R = 1;
  dH[t_] := -2 \[CapitalDelta]*Sin[2*t]; 
  ygrid = Range[-1, 1, dy]; (* Updated grid to include -1 to 1 *)
  ny = Length[ygrid];
  
  dy2 = NDSolve`FiniteDifferenceDerivative[Derivative[2], ygrid, 
     DifferenceOrder -> 4]@"DifferentiationMatrix";
  dy1 = NDSolve`FiniteDifferenceDerivative[Derivative[1], ygrid, 
     DifferenceOrder -> 4]@"DifferentiationMatrix"; 
  fvar0 = ConstantArray[0, {ny, nmax}]; (* Initial condition array for V *)
  
  Do[
   (* Extract the solution for the current time step *)
   fvar = fvar0[[All, i]];
   
   (* Define the operators based on 2.5 and 2.6 *)
   LV = (1/H[(i + 1) dt])*(2 H[(i + 1) dt] dy2.fvar + 
        ygrid H[(i + 1) dt] (dy1.fvar) + (1/(R H[(i + 1) dt]))*fvar);
   
   NV = (1/H[(i + 1) dt])*(fvar (dy1.fvar));
   
   (* The main equation V_eta_eta_t = LV + NV *)
   eqf = (fvar0[[All, i + 1]] - fvar)/dt - (LV + NV);
   
   (* Boundary Conditions based on Eq 2.4 *)
   eqf[[1]] = fvar[[1]] - (-dH[(i + 1) dt]); (* V(-1, t) = -dH *)
   eqf[[-1]] = fvar[[-1]] - dH[(i + 1) dt];  (* V(1, t) = dH *)
   
   eqf[[1]] = (dy1.fvar)[[1]];               (* V_eta(-1, t) = 0 *)
   eqf[[-1]] = (dy1.fvar)[[-1]];             (* V_eta(1, t) = 0 *)
   
   eqns = eqf;
   {vec, mat} = CoefficientArrays[eqns, fvar];
   
   (* Solve the linear system *)
   sol = LinearSolve[mat, -vec]; 
   
   (* Update the solution array for the next step *)
   fvar0[[All, i + 1]] = sol;, 
   {i, nmax - 1}]; 
  
  fvar0
]

delta = 0.5; (* Example delta value; you can choose other values too *)
result = soln[delta];

(* Plot the solution at the last time step *)
ListLinePlot[
 Transpose[{ygrid, result[[All, -1]]}], 
 Frame -> True, 
 FrameLabel -> {"y", "V(y, t_final)"}, 
 PlotLabel -> "Final Solution for V(y, t) at t = nmax * dt"
]

Range for $\Delta$ is: $0 < \Delta < 1$

They have used the Crank-Nicolson method.

The Manuscript with title "The Onset of Chaos in a class of Navier-Stokes Solutions" is also provided here:

https://doi.org/10.1017/S0022112099005364

The above stated PDE is on page number $61$ of the given manuscript against the equation number $(2.4)$ along with linear and nonlinear terms presented in equations $(2.5)$ and $(2.6)$ respectively. The code is complete in all aspect but still not working.

Need your guideline.

If you need more information or think I have missed something, please let me know.