How to fixed NDSolve error by Mathematica Program

Question

This is non-linaer system of wave equation. To solve this equation I used NDSolve (Mathematica). I try to solved non-linear wave equation numerically by using FDM and NDSolve but I got some error. I think the code before solf=NDSolve is ok but may be from solf is creat problem, I could not get solution. In this program I used multiple amplitude and frequency, the program code and which error sappear are given below:

(Errors) enter image description here

`

Clear[t, \[Delta], R, An, Q0, A, v, \[Nu], dt, Nf, cn, \[Omega], \
Eq1D, Eq2D, Dt1, Dt2, Dt3, Dt4, Tp, f, n, nn, s, xf, LBC, LDEq, Lv, \
LA, Lvsol, LAsol, solN]

QAmplist = {1.2879968521168974`*^-7 - 3.032561883723498`*^-8 I, 
   1.3196528666022895`*^-7 - 5.591765158902973`*^-8 I, 
   3.69285439379095`*^-7 - 5.822321853105695`*^-7 I, 
   1.5819124403198096`*^-8 + 3.1965417764762`*^-7 I};
fplist = {1/5, 2/5, 3/5, 4/5};

AmpQ1 = {};
AmpQ2 = {};
TU = {};
TV = {};
TP = {};
TQ = {};
Print["Start"]
Table[
  cn = 111 ;(*m/s*)
  R = 4 10^-3; (*m*)
  Q0 = 2.03 10^-3/60; (*m^3/s*)
  \[Mu] = 0.0160; (*pa-s*)
  \[Rho] = 1167;(*kg/m^3*)
  \[Nu] = \[Mu] /\[Rho];
  L = 10;
  \[Delta] = 0.1;
  An = \[Pi] R^2;
  (*Tp=0.2;\[Omega]=2 \[Pi] 1/Tp;*)
  Tp = 1/fplist[[f]];
  \[Omega] = 2 \[Pi] 1/Tp;
  (*Tp=0.2;\[Omega]=2 \[Pi] 1/Tp;*)
  
  (*Print["Disharge fluctuation=",\[Delta] Q0];*)
  
  (* Assuming a time-periodic solution and solving with the finite-
  difference method*)
  
  (* boundary condition：QI[s_]=Q0(1+\[Delta] Sin[\[Omega] s]); *)
  (* the period is Tp/2 *)
  (*(*QI[t_]=QAmplist(*)*( (Abs[Sin[\[Omega] s/2]]+Abs[Sin[\[Omega] s/
  2+\[Pi]/2]]));**)
  (*Print@Plot[QI[s],{s,0,\[Omega]},PlotRange\[Rule]{0, 
  1.5Q0 QAmplist},Frame\[Rule]True,FrameLabel\[Rule]{"t(s)",
  "Q(m^3/s)"}];*)
  
  Nf = 30; dt = Tp/Nf;
  t[n_] = (n - 1) Tp/Nf/2;
  
  Dt1[f_] := (f[n + 1][x] - f[n - 1][x])/(2 dt);
  (*Dt2[f_]:=(f[n+1][x]-2f[n][x]+f[n-1][x])/(dt^2);
  Dt3[f_]:=(f[n+2][x]-2f[n+1][x]+2f[n-1][x]-f[n-2][x])/(2dt^3); 
  Dt4[f_]:=(f[n+2][x]-4f[n+1][x]+6f[n][x]-4f[n-1][x]+f[n-2][
  x])/(dt^4); *)
  (* LDt5[f_]:=(-f[-3 +n][x]+4 f[-2 +n][x]-5 f[-1+n][x]+5 f[1+n][
  x]-4 f[2 +n][x]+f[3 +n][x])/(2dt^5); *)
  
  (* discretized as a time period function: n=1,...,Nf *)
  Eq1D[n_] = Dt1[A] + v[n][x] D[A[n][x], x] + A[n][x] D[v[n][x], x];
  Eq2D[n_] = 
   4/3 A[n][x] Dt1[v] + (An^2 cn^2 D[A[n][x], x])/A[n][x]^2 + 
    v[n][x] (8 \[Pi] \[Nu] + A[n][x] D[v[n][x], x]);
  (* time periodic condition*)
  Do[
   A[1 - nn][x_] = A[Nf + 1 - nn][x];
   v[1 - nn][x_] = v[Nf + 1 - nn][x];
   A[Nf + nn][x_] = A[nn][x];
   v[Nf + nn][x_] = v[nn][x];
   , {nn, 1, 2}];
  
  (* boundary condition *)
  
  LBC = Flatten[ 
    Table[{A[n][0] v[n][0] == 
       QAmplist [[f]] ( 1 + Abs[Sin[\[Omega] t[n]]]) (*QI[t[n]]*), 
      A[n][L] == An}, {n, 1, Nf}] ];
  LDEq = Flatten[ Table[{Eq1D[n] == 0, Eq2D[n] == 0}, {n, 1, Nf}] ];
  Lv[x_] = Table[v[n][x], {n, 1, Nf}];
  LA[x_] = Table[A[n][x], {n, 1, Nf}];
  LQ[x_] = Table[A[n][x] v[n][x], {n, 1, Nf}];
  
  solf = 
   NDSolve[Flatten[{LDEq, LBC}], Flatten[{Lv[x], LA[x]}], {x, 0, 5}];
  
  , {f, 1, Length[fplist]}];
Print["finish"];

`