This Mathematica notebook demonstrates the numerical method described in the Appendix of the article  P. Bizoń & P. Mach, `GLOBAL DYNAMICS OF A YANG-MILLS FIELD ON AN ASYMPTOTICALLY HYPERBOLIC SPACE’

We solve the initial value problem for the equation
(∂_(ττ))W + 2x (∂_(τx))W + (∂_τ)W=(∂_x)((1 - (x^2))(∂_x)W)+ℓ(ℓ+1)W(1-(W^2)).

Set the number of Chebyshev polynomials = nmax+1

nmax=30

Set the value of ℓ

L=3+1/2

Construct the equations for expansion coefficients a[n][t]

Define the cutoff function for the coefficients

cutoff[n_]:=If[n<0,0,If[n>nmax,0,1]]

cuta[n_]:=cutoff[n]a[n][t]

Compute the coefficients for the (W^3)term

cubic[k_]:=If[k>0,Sum[(1/4)(cuta[k-m-n]+cuta[k+m+n]+cuta[-k+m+n]+cuta[k-Abs[m-n]]+cuta[k+Abs[m-n]]+cuta[-k+Abs[m-n]])a[m][t]a[n][t],{m,0,nmax},{n,0,nmax}],(1/4)(a[0][t])^3+Sum[(1/4)(cuta[m+n]+cuta[Abs[m-n]]+cuta[-Abs[m-n]])a[m][t]a[n][t],{m,0,nmax},{n,0,nmax}]]

Compute the coefficients for the x (∂_x)W term

z[k_]:=If[k==0,Sum[n a[n][t],{n,0,nmax,2}],k a[k][t]+2Sum[n a[n][t],{n,k+2,nmax,2}]]

Compute the coefficients for the x (∂_(τx))W term

y[k_]:=If[k==0,Sum[n a[n]'[t],{n,0,nmax,2}],k a[k]'[t]+2Sum[n a[n]'[t],{n,k+2,nmax,2}]]

Define the set of equations for the expansion coefficients

eqns=Table[a[n]''[t]+a[n]'[t]+(n^2 - L(L+1))a[n][t]+2y[n]+z[n]+L(L+1) cubic[n]==0,{n,0,nmax}];

Define the set of variables (unknowns)

vars=Table[a[n][t],{n,0,nmax}];

Define initial conditions

ic1a={a[0][0]==0,a[1][0]==2,a[2][0]==0,a[3][0]==0,a[4][0]==0,a[5][0]==1,a[6][0]==0}

ic1b=Table[a[n][0]==0,{n,7,nmax}];

Define initial conditions for the momentum

ic2=Table[a[n]'[0]==0,{n,0,nmax}];

Collect all equations together

eqnstot=Join[eqns,ic1a,ic1b,ic2];

Solve the equations. Note that in some cases one should work using much higher computational precision (WorkingPrecision → 30 or so)

sol=NDSolve[eqnstot,vars,{t,0,50}];

Plot the time dependence of the coefficients

Plot[Evaluate[vars/.sol],{t,0,50},PlotRange->{{0,50},{-2,2}},AxesLabel->{t,(Subscript[a, n])[t]}]

Construct the solution W

ww=Sum[a[n][t]ChebyshevT[n,x],{n,0,nmax}]/.sol;

Plot the solution

Animate[Plot[ww/.t->T,{x,-1,1},PlotRange->{{-1,1},{-3,3}},ImageSize->700,AxesLabel->{x,W}],{T,0,20}]