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={}},Δ=delta;H[t_]:=1+Δ*Cos[2*t];dH[t_]:=-2*Δ*Sin[2*t];(*R=100;*)(*Createthegridanddifferentiationmatrices*)ygrid=Range[0,1,dy];ny=Length[ygrid];dy2=NDSolve`FiniteDifferenceDerivative[Derivative[2],ygrid,DifferenceOrder4]@"DifferentiationMatrix";dy1=NDSolve`FiniteDifferenceDerivative[Derivative[1],ygrid,DifferenceOrder4]@"DifferentiationMatrix";(*InitializevariablesforFandG*){fvar1,gvar1}=MakeVariables[{ff,gg},ny];fvar0=ConstantArray[0,{ny,nmax}];gvar0=ConstantArray[0,{ny,nmax}];(*Time-steppingloop*)Do[(*UpdatevariablesforFandGatthecurrenttimestep*){fvar,gvar}={fvar0[[All,i]],gvar0[[All,i]]};(*Definetheequations*)eqf=dy2.fvar1+(dy1.fvar1)/ygrid-fvar1/ygrid^2+gvar1H[(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/(ygridH[(i+1)dt])fvar*gvar1-(1/(RH[(i+1)dt]^2))*((dy2.gvar1)+(dy1.gvar1)/ygrid-gvar1/ygrid^2)//Quiet;(*Applyboundaryconditions*)eqf[[1]]=fvar1[[1]];(*ensuresthatF(0,t)=0aty=0*)eqf[[-1]]=fvar1[[-1]]+dH[(i+1)dt];(*ensuresthatF(1,t)=-H'[t]aty=1*)eqg[[1]]=gvar1[[1]];(*ensuresthatG(0,t)=0att=0*)eqg[[-1]]=(dy1.fvar1)[[-1]]-dH[(i+1)dt];(*ensuresthatF'(1,t)=H'[t]aty=1*)(*Solvethesystem*)eqns=Join[eqf,eqg];var=Join[fvar1,gvar1];{vec,mat}=CoefficientArrays[eqns,var];sol=LinearSolve[mat,-vec];(*ExtractsolutionsforFandGatthenexttimestep*){fvar0[[All,i+1]],gvar0[[All,i+1]]}=Partition[sol,ny];(*InterpolateG(y,t)andstoreG(1/4,t)*)gAtTarget=Interpolation[Transpose[{ygrid,gvar0[[All,i+1]]}]];(*StoreG(1/4,t)onlyfortfrom8100to8200*)If[i≥8100&&i≤8201,AppendTo[gValues,{i,gAtTarget[yTarget]}];];,{i,1,nmax-1}];(*ReturnthesolutionsforF,GandthestoredGvalues*){fvar0,gvar0,gValues}];(*Callthefunctionfordelta=0.2andstoreG(1/4,t)from8100to8200*){fSol,gSol,gValues}=soln[.2,2050];//AbsoluteTiming;gValues;(*PlotthevaluesofG(1/4,t)*)sajjad=ListLinePlot[gValues,GridLinesAutomatic,ImageSizeLarge,FrameAllTrue,FrameLabel{"t","G(1/4,t)"},PlotStyle{Red,Thick}]
