Hello everyone. I am new in MATHEMATICA. Would anyone help/guide how to solve the following system of coupled PDE’s with time.? I have tried to solve but did not get the desire results (there is no smoothness in the graph). Firstly, I want to draw the fig 10 a in paper. Then fig 8 etc. The code has been attached. The following system got from this paper with DOI, (DOI: 10.1017/S0022112003003835) Thank you.

G_t-H^'/H*η*G_η-1/H*F*G_η+1/H*F_η*G+(2*F*G)/(η*H)=1/(H^2*R) [G_ηη+G_η/η-G/η^2 ]

F_ηη+F_η/η-F/η^2 =-H^2*G
With:
                 Initial conditions
F(η,0)=0,                 G(η,0)=0
                 Boundary conditions
F(0,t)=0                       G(0,t)=0,
F(1,t)=-H'                F_η (1,0)=H'

H(t)=1+Δ*cos⁡(2t)
Where               0<Δ<1

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

soln[delta_, Renolds_] := 
  Module[{fvar0, gvar0, fvar1, gvar1, dy1, dy2, ny, ygrid, 
    R = Renolds, dt = 1, nmax = 8201, dy = 1/10, yTarget = 1/4, 
    gAtTarget, gValues = {}}, \[CapitalDelta] = delta;
   H[t_] := 1 + \[CapitalDelta]*Cos[2*t];
   dH[t_] := -2*\[CapitalDelta]*Sin[2*t];
   (*R=100;*)
   (*Create the grid and differentiation matrices*)
   ygrid = Range[0, 1, dy];
   ny = Length[ygrid];
   dy2 = NDSolve`FiniteDifferenceDerivative[Derivative[2], ygrid, 
      DifferenceOrder -> 4]@"DifferentiationMatrix";
   dy1 = NDSolve`FiniteDifferenceDerivative[Derivative[1], ygrid, 
      DifferenceOrder -> 4]@"DifferentiationMatrix";
   (*Initialize variables for F and G*){fvar1, gvar1} = 
    MakeVariables[{ff, gg}, ny];
   fvar0 = ConstantArray[0, {ny, nmax}];
   gvar0 = ConstantArray[0, {ny, nmax}];
   (*Time-stepping loop*)
   Do[(*Update variables for F and G at the current time step*){fvar, 
      gvar} = {fvar0[[All, i]], gvar0[[All, i]]};
    (*Define the equations*)
    eqf = dy2.fvar1 + (dy1.fvar1)/ygrid - fvar1/ygrid^2 + 
       gvar1 H[(i + 1) dt]^2 // Quiet;
    eqg = (gvar1 - gvar)/dt - 
       dH[(i + 1) dt]/
         H[(i + 1) dt] ygrid (dy1.gvar1) - (dy1.gvar1) fvar/
         H[(i + 1) dt] + 1/H[(i + 1) dt] (dy1.fvar) gvar1 + 
       2/(ygrid H[(i + 1) dt]) fvar*
        gvar1 - (1/(R H[(i + 1) dt]^2))*((dy2.gvar1) + (dy1.gvar1)/
           ygrid - gvar1/ygrid^2) // Quiet;
    (*Apply boundary conditions*)
    eqf[[1]] = fvar1[[1]];(*ensures that F(0,t)=0 at y=0*)
    eqf[[-1]] = fvar1[[-1]] + dH[(i + 1) dt];(*ensures that F(1,
    t)=-H'[t] at y=1*)eqg[[1]] = gvar1[[1]];(*ensures that G(0,t)=
    0 at t=0*)
    eqg[[-1]] = (dy1.fvar1)[[-1]] - 
      dH[(i + 1) dt];(*ensures that F'(1,t)=H'[t] at y=
    1*)(*Solve the system*)eqns = Join[eqf, eqg];
    var = Join[fvar1, gvar1];
    {vec, mat} = CoefficientArrays[eqns, var];
    sol = LinearSolve[mat, -vec];
    (*Extract solutions for F and G at the next time step*)
    {fvar0[[All, i + 1]], gvar0[[All, i + 1]]} = Partition[sol, ny];
    (*Interpolate G(y,t) and store G(1/4,t)*)
    gAtTarget = Interpolation[Transpose[{ygrid, gvar0[[All, i + 1]]}]];
    (*Store G(1/4,t) only for t from 8100 to 8200*)
    If[i >= 8100 && i <= 8201, 
     AppendTo[gValues, {i, gAtTarget[yTarget]}];];, {i, 1, nmax - 1}];
   (*Return the solutions for F,G and the stored G values*)
   {fvar0, gvar0, gValues}];

(*Call the function for delta=0.2 and store G(1/4,t) from 8100 to \
8200*)
{fSol, gSol, gValues} = soln[.2, 2050]; // AbsoluteTiming;

gValues;

(*Plot the values of G(1/4,t)*)


sajjad = ListLinePlot[gValues, GridLines -> Automatic, 
  ImageSize -> Large, Frame -> AllTrue, 
  FrameLabel -> {"t", "G(1/4,t)"}, PlotStyle -> {Red, Thick}]

Attachments:

Fig 10 a.nb