Método de diferencias finitas para la Ecuación de Boussinesq mejorada

La Ecuación de Boussinesq mejorada es: (Subscript[u, tt])-(Subscript[u, xx])-(Subscript[u, xxtt])-(Subscript[((u^2)), xx])= 0.

({{ (Subscript[u, tt])-(Subscript[u, xx])-(Subscript[u, xxtt])-(Subscript[((u^2)), xx])= 0,             for      0⩽x⩽L,    t≥0,
u(x,0)=f(x),                                                  for     a⩽x⩽b,  
(Subscript[u, t])(x,0)=g(x),                                               for      a⩽x⩽b,}, {u(a,t)= u(b,t)=0,                                    for     t>0
(Subscript[u, x]) (a,t)=(Subscript[u, x])(b,t)=0,                             for     t>0
(Subscript[u, xx]) (a,t)=(Subscript[u, xx])(b,t)=0,                         for     t>0}})

PRIMER EJEMPLO : UN SOLITÓN

Clear["Global`*"]

(*Definimos los parámetros*)

A=0.5;
c=(Sqrt[1+((2*A)/3)]);
x0=80;

(*Definimos u(x,t) y sus derivadas*)

u[x_,t_]:=u[x,t]=A*(((Sech[(1/c)*(Sqrt[(A/6)])*(x-c*t-x0)]))^2);
ux[x_,t_]:=ux[x,t]=D[u[x,t],x];
uxx[x_,t_]:=uxx[x,t]=D[u[x,t],{x,2}]

(*Dibujamos la gráfica de u(x,t) para ver si se cumplen las condiciones de contorno*)

Plot3D[u[x,t],{x,0,200},{t,0,70},PlotStyle->Cyan,Mesh->None,AxesLabel->Automatic, PlotRange->Full,PlotPoints->200,PlotLabel->"u(x,t)"]

Plot[u[x,0],{x,0,200},PlotStyle->Cyan,Mesh->None,AxesLabel->Automatic, PlotRange->Full,PlotLabel->"u(x,0)"]

Plot[u[x,70],{x,0,200},PlotStyle->Cyan,Mesh->None,AxesLabel->Automatic, PlotRange->Full,PlotLabel->"u(x,70)"]

(*Comprobamos analíticamente que se cumple que las 3 condiciones de cortono (CC) son cero*)

u[0,70]
u[200,70]

ux[x,t]/.{x->0,t->70}
ux[x,t]/.{x->200,t->70}

uxx[x,t]/.{x->0,t->70}
uxx[x,t]/.{x->200,t->70}

(*Hacemos el programa que resuelve el problema de la ecuación mejorada de Boussinesq (3.1) para el primer ejemplo*)

Clear["Global`*"];
t0=TimeUsed[ ];

L=200; 
M=1999; (*19999 puntos interiores -> 20000 intervalos*)
A=0.5; (*Amplitud de la onda*)
c=(Sqrt[1+((2*A)/3)]); (*Velocidad de la onda*)
x0=80;
Δx = L/(M+1); (*El paso del espacio: k*)
tfin=70.0;
Δt=0.1; (*El paso del tiempo: h*)
niteraciones=Ceiling[tfin/Δt,1]; 
(*Número de iteraciones en tiempo. Ceiling da el menor múltiplo de 1 mayor o igual a tfin/Δt*)
α1=1;α2=1;α3=-1; (*Valores de alfa*)

Print["\!\(TraditionalForm\`\*SubscriptBox[\(u\), \(tt\)] - \*SubscriptBox[\(u\), \(xx\)] - \*SubscriptBox[\(u\), \(xxtt\)] - \*SubscriptBox[\((\*SuperscriptBox[\(u\), \(2\)])\), \(xx\)] = \\ 0\)",",   A = ",  A,",   c = ",  c,",   M = ", M, ",  Δx = ",  Δx,",   Δt = ",  Δt, ",   niteraciones = ",niteraciones, ",   real final time = ", N[niteraciones Δt]]

(*Definimos las condiciones iniciales (CI)*)
f[x_]:=f[x]=A*(((Sech[(1/c)*(Sqrt[(A/6)])*(x-x0)]))^2);(*u(x,0)*)

g1=Plot[f[x],{x,0,L},PlotRange-> All,PlotLegends->{"Original function f(x)"},PlotStyle-> Cyan]

u[x_,t_]:=u[x,t]=A*(((Sech[(1/c)*(Sqrt[(A/6)])*(x-c*t-x0)]))^2);
(*Aproximamos la primera CI*)
U[s_,0]:= U[s,0] =f[x]/. x -> s Δx ; 

Ut[s_,0]:=Ut[s,0]=D[u[x,t],t]/. x -> s Δx /. t -> 0
Utt[s_,0]:=Utt[s,0]=D[u[x,t],{t,2}]/. x -> s Δx /. t -> 0
(*Aproximamos la segunda CI*) 
U[s_,1]:= U[s,1] =N[U[s,0]+Δt*Ut[s,0]+((((Δt))^2)/2)*Utt[s,0]];

(*Condiciones de contorno*)
U[0,r_]:=0; U[M+1,r_]:=0;
Ux[0,r_]:=0;Ux[M+1,r_]:=0;
Uxx[0,r_]:=0;Uxx[M+1,r_]:=0;

U[-1,r_]:=0; U[M+2,r_]:=0;
Ux[-1,r_]:=0;Ux[M+2,r_]:=0;
Uxx[-1,r_]:=0;Uxx[M+2,r_]:=0;

(*Definimos 3.3 y 3.4*)
Ux[s_,0]:=Ux[s,0]=(1/(12*Δt))(-U[s+2,0]+8*U[s+1,0]-8*U[s-1,0]+U[s-2,0]);
Uxx[s_,0]:=Uxx[s,0]=(1/(12*(((Δt))^2)))(-U[s+2,0]+16*U[s+1,0]-30*U[s,0]+16*U[s-1,0]-U[s-2,0]);
Ux[s_,1]:=Ux[s,1]=(1/(12*Δt))(-U[s+2,1]+8*U[s+1,1]-8*U[s-1,1]+U[s-2,1]);
Uxx[s_,1]:=Uxx[s,1]=(1/(12*(((Δt))^2)))(-U[s+2,1]+16*U[s+1,1]-30*U[s,1]+16*U[s-1,1]-U[s-2,1]);

(*Definimos los lambdas del sistema 3.6*)
l1[s_,r_]:=l1[s,r]=1-2*((α1)^2)*(((Δx))^2)*Uxx[s,r]-2*α1*α2*(((Δx))^2)*Uxx[s,r-1]-2*α1*α3*(((Δx))^2)*Uxx[s,r-2];
l2[s_,r_]:=l2[s,r]=-2*((α1)^2)*(((Δx))^2)*Ux[s,r]-4*α1*α2*(((Δx))^2)*Ux[s,r-1]-4*α1*α3*(((Δx))^2)*Ux[s,r-2];
l3[s_,r_]:=l3[s,r]=-1-α1*(((Δx))^2)-2*α1*α2*(((Δx))^2)*U[s,r-1]-2*α1*α3*(((Δx))^2)*U[s,r-2];
l4[s_,r_]:=l4[s,r]=(2+2*((α2)^2)*(((Δx))^2)*Uxx[s,r-1]+2*α2*α3*(((Δx))^2)*Uxx[s,r-2])*U[s,r-1]+
(-1+2*α2*α3*(((Δx))^2)*Uxx[s,r-1]+2*((α3)^2)*(((Δx))^2)*Uxx[s,r-2])*U[s,r-2]+
(2*((α2)^2)*(((Δx))^2)*Ux[s,r-1]+4*α2*α3*(((Δx))^2)*Ux[s,r-2])*Ux[s,r-1]+
2*((α3)^2)*(((Δx))^2)*Ux[s,r-2]*Ux[s,r-2]+(α2*(((Δx))^2)-2)*Uxx[s,r-1]+(α3*(((Δx))^2)+1)*Uxx[s,r-2];

W4[r_]:=W4[r]=Table[l4[i,r],{i,1,M}];

(*Definimos la matriz para resolver el sistema 3.6*)

T[r_]:=T[r]=SparseArray[{i_,j_}/;i==j+2->(((l2[i,r])/(12*(Δt)))-((l3[i,r])/(12*(((Δt))^2)))),{M,M}]+
SparseArray[{i_,j_}/;i==j+1->(-((2*l2[i,r])/(3*(Δt)))+((4*l3[i,r])/(3*(((Δt))^2)))),{M,M}]+
SparseArray[{i_,j_}/;i==j->(l1[i,r]-((5*l3[i,r])/(2*(((Δt))^2)))),{M,M}]+
SparseArray[{i_,j_}/;i==j-1->(((2*l2[i,r])/(3*(Δt)))+((4*l3[i,r])/(3*(((Δt))^2)))),{M,M}]+
SparseArray[{i_,j_}/;i==j-2->(-((l2[i,r])/(12*(Δt)))-((l3[i,r])/(12*(((Δt))^2)))),{M,M}];

(*Proceso de refinamiento iterativo para obtener (ux)r+1,s y (uxx)r+1,s*)
refinamiento[N0_,TOL_,M_, Ux0_,Ux1_,AuxUx2_, Uxx0_,Uxx1_,AuxUxx2_, α1_,α2_,α3_, iter_]:=
Module[{N0var= N0, Tolvar=TOL,Mvar=M,Ux0var = Ux0,Ux1var = Ux1,AuxUx2var = AuxUx2,Uxx0var = Uxx0,Uxx1var = Uxx1,AuxUxx2var = AuxUxx2,α1var=α1,α2var =α2,α3var=α3, itervar=iter,i,s,XX,X,Δx, Sol, UU, j,AuxUUx2var,AuxUUxx2var,AuxMat, AuxVect1, AuxVect2, AuxVect,l1,l2,l1new,l2new},
Δx = 200/(Mvar+1);
i=1;
XX=Table[X[s],{s,1, 2 M}];
l1=1-2*((α1)^2)*(((Δx))^2)*Uxx1var-2*α1*α2*(((Δx))^2)*Uxx1var-2*α1*α3*(((Δx))^2)*Uxx0var;
l2=-2*((α1)^2)*(((Δx))^2)*Ux1var-4*α1*α2*(((Δx))^2)*Ux1var-4*α1*α3*(((Δx))^2)*Ux0var;

AuxMat := AuxMat= SparseArray[{i_,j_}/;i==j->(-2*((α1var)^2)*(((Δx))^2)),{2 M,2 M}];
AuxVect1:= AuxVect1= 1-2*α1var*α2var*(((Δx))^2)*Uxx1var-2*α1var*α3var*(((Δx))^2)*Uxx0var-l1;

AuxVect2:= AuxVect2= -4*α1var*α2var*(((Δx))^2)*Ux1var-4*α1var*α3var*(((Δx))^2)*Ux0var - l2;

AuxVect=Join[AuxVect1, AuxVect2];

While[i≤N0var,
itervar=i;
Print["itervar = ",itervar];
Sol=NSolve[{AuxMat.XX==AuxVect},XX];
UU= XX/. Sol;


For[j=1,j≤ M,j++, AuxUUxx2var[j]= UU[[1,j]]];
For[j=1,j≤ M,j++, AuxUUx2var[j]= UU[[1,j+M]]];
AuxUx2var=Table[AuxUUx2var[s],{s,1,M}];
AuxUxx2var=Table[AuxUUxx2var[s],{s,1,M}];

If[Norm[AuxUxx2var-Uxx1var]<Tolvar && Norm[AuxUx2var-Ux1var]<Tolvar ,
{Print[" Control  YES"];
Return[AuxUx2var,Module];
Return[AuxUxx2var,Module];
Return[itervar ,Module]};
];

Print[" Control  NO"];

l1new=1-2*((α1)^2)*(((Δx))^2)*AuxUxx2var-2*α1*α2*(((Δx))^2)*Uxx1var-2*α1*α3*(((Δx))^2)*Uxx0var;
l2new=-2*((α1)^2)*(((Δx))^2)*AuxUx2var-4*α1*α2*(((Δx))^2)*Ux1var-4*α1*α3*(((Δx))^2)*Ux0var;


AuxVect1:= AuxVect1= 1-2*α1var*α2var*(((Δx))^2)*Uxx1var-2*α1var*α3var*(((Δx))^2)*Uxx0var-l1new;

AuxVect2:= AuxVect2= -4*α1var*α2var*(((Δx))^2)*Ux1var-4*α1var*α3var*(((Δx))^2)*Ux0var - l2new;

AuxVect=Join[AuxVect1, AuxVect2];

i++];

Return[AuxUx2var,Module];
Return[AuxUxx2var,Module];
Return[itervar ,Module];

]

N0=50;
TOL=0.0000001;
XX=Table[X[s],{s,1, M}];

(*Hacemos dos bucles para resolver nuestro sistema para cada valor de s y cada valor de r*)
For[r=2,r≤niteraciones,r++,     (*AQUÍ   EMPIEZA EL BUCLE TEMPORAL*)
{Print[" r = ", r],
iter=0;

For[s=1, s≤ M,s++,{
UU0[s]= U[s,r-2];
UUx0[s]= Ux[s,r-2];
UUxx0[s]= Uxx[s,r-2];

UU1[s]= U[s,r-1];
UUx1[s]= Ux[s,r-1];
UUxx1[s]= Uxx[s,r-1];

UU2[s]= U[s,r-1];
UUx2[s]= Ux[s,r-1];
UUxx2[s]= Uxx[s,r-1];
}];

UUx2[-1]:=0;UUx2[M+2]:=0;
UUxx2[-1]:=0;UUxx2[M+2]:=0;

UUUx0=Table[Ux[i,r-2],{i,1,M}];
UUUxx0=Table[Uxx[i,r-2],{i,1,M}];
UUUx1=Table[Ux[i,r-1],{i,1,M}];
UUUxx1=Table[Uxx[i,r-1],{i,1,M}];
UUUx2=Table[Ux[i,r-1],{i,1,M}];
UUUxx2=Table[Uxx[i,r-1],{i,1,M}];

refinamiento[N0,TOL,M,UUUx0,UUUx1,UUUx2, UUUxx0,UUUxx1,UUUxx2,α1,α2,α3, iter];


For[j=1,j≤ M,j++, Ux[j,r]=UUUx2[[j]]];
For[j=1,j≤ M,j++, Uxx[j,r]=UUUxx2[[j]]];


W4[r]=Table[l4[i,r],{i,1,M}];
UU:=UU=Table[0,{i,1,M}];
AA= T[r];
bb= W4[r];

XX=Table[X[s],{s,1, M}];
Sol=Solve[AA.XX== bb,XX];
UU= XX/. Sol;


For[s=1,s≤ M,s++, U[s,r]= UU[[1,s]]];

For[s=1,s≤ M,s++, Ux[s,r]=(1/(12*Δt))(-U[s+2,r]+8*U[s+1,r]-8*U[s-1,r]+U[s-2,r])];
For[s=1,s≤ M,s++, Uxx[s,r]=(1/(12*(((Δt))^2)))(-U[s+2,r]+16*U[s+1,r]-30*U[s,r]+16*U[s-1,r]-U[s-2,r])];

}
 ];    (*AQUÍ TERMINA EL BUCLE TEMPORAL*)

Plot1=Table[ListPlot[Table[{Δx s, U[s,r]},{s,0,M}], Joined-> True, PlotRange-> All,PlotStyle->Cyan, AxesLabel -> {"x", " "}, PlotLabel -> "T[x,t], t=" <> ToString[Δt r]], {r,0,niteraciones,25}]

Data3D = Table[{Δx s,0 , U[s,0]},{s,0,M+1}];
Do[Data3D=Join[Data3D,Table[{Δx s,Δt r , U[s,r]},{s,0,M+1}]],{r,1,niteraciones,10}]; 
Plot2=ListPlot3D[Data3D,PlotRange-> All,Mesh->{0,0},AxesLabel->{"Space - x","Time - t","u"},PlotLabel ->"Boussinesq improved equation",PlotStyle->Cyan]

t1=TimeUsed[ ];
Print[" El tiempo total empleado son ", t1-t0, " segundos."]

SEGUNDO EJEMPLO : INTERACCIÓN ENTRE DOS SOLITONES

Clear["Global`*"]

(*Definimos los parámetros*)

A1=0.4; A2=0.2;
c1=(Sqrt[1+((2*A1)/3)]);c2=-(Sqrt[1+((2*A2)/3)]);
x01=90;x02=130;

(*Definimos u(x,t) y sus derivadas*)

u[x_,t_]:=u[x,t]=A1*(Sech[(1/c1)*(Sqrt[(A1/6)])*(x-c1*t-x01)])+A2*(Sech[(1/c2)*(Sqrt[(A2/6)])*(x-c2*t-x02)]);
ux[x_,t_]:=ux[x,t]=D[u[x,t],x];
uxx[x_,t_]:=uxx[x,t]=D[u[x,t],{x,2}];

(*Dibujamos la gráfica de u(x,t) para ver si se cumplen las condiciones de contorno*)

Plot3D[u[x,t],{x,0,220},{t,0,60},PlotStyle->Cyan,Mesh->None,AxesLabel->Automatic, PlotRange->Full,PlotPoints->200,PlotLabel->"u(x,t)"]

Plot[u[x,0],{x,0,220},PlotStyle->Cyan,Mesh->None,AxesLabel->Automatic, PlotRange->Full,PlotLabel->"u(x,0)"]

Plot[u[x,60],{x,0,220},PlotStyle->Cyan,Mesh->None,AxesLabel->Automatic, PlotRange->Full,PlotLabel->"u(x,60)"]

(*Comprobamos que se cumplen las 3 CC*)

u[0,60]
u[220,60]

ux[x,t]/.{x->0,t->60}
ux[x,t]/.{x->200,t->60}

uxx[x,t]/.{x->0,t->60}
uxx[x,t]/.{x->200,t->60}

(*Hacemos el programa que resuelve el problema de la ecuación mejorada de Boussinesq (3.1) para el primer ejemplo*)

Clear["Global`*"];
t0=TimeUsed[ ];

L=200; 
M=1999; (*19999 puntos interiores -> 20000 intervalos*)
A1=0.4; A2=0.2; (*Amplitud de la onda*)
c1=(Sqrt[1+((2*A1)/3)]);c2=-(Sqrt[1+((2*A2)/3)]); (*Velocidad de la onda*)
x01=90;x02=130;
Δx = L/(M+1); (*El paso del espacio: k*)
tfin=60.0;
Δt=0.1; (*El paso del tiempo: h*)
niteraciones=Ceiling[tfin/Δt,1]; 
(*Número de iteraciones en tiempo. Ceiling da el menor múltiplo de 1 mayor o igual a tfin/Δt*)
α1=1;α2=1;α3=-1; (*Valores de alfa*)

Print["\!\(TraditionalForm\`\*SubscriptBox[\(u\), \(tt\)] - \*SubscriptBox[\(u\), \(xx\)] - \*SubscriptBox[\(u\), \(xxtt\)] - \*SubscriptBox[\((\*SuperscriptBox[\(u\), \(2\)])\), \(xx\)] = \\ 0\)",",   A = ",  A,",   c = ",  c,",   M = ", M, ",  Δx = ",  Δx,",   Δt = ",  Δt, ",   niteraciones = ",niteraciones, ",   real final time = ", N[niteraciones Δt]]

(*Definimos las condiciones iniciales (CI)*)
u[x_,t_]:=u[x,t]=A1*(Sech[(1/c1)*(Sqrt[(A1/6)])*(x-c1*t-x01)])+A2*(Sech[(1/c2)*(Sqrt[(A2/6)])*(x-c2*t-x02)]);
ut[x_,t_]:=ut[x,t]=D[u[x,t],t];
f[x_]:=f[x]=u[x,0];(*u(x,0)*)

g1=Plot[f[x],{x,0,L},PlotRange-> All,PlotLegends->{"Original function f(x)"},PlotStyle-> Cyan]

(*Aproximamos la primera CI*)
U[s_,0]:= U[s,0] =f[x]/. x -> s Δx ; 

Ut[s_,0]:=Ut[s,0]=D[u[x,t],t]/. x -> s Δx /. t -> 0
Utt[s_,0]:=Utt[s,0]=D[u[x,t],{t,2}]/. x -> s Δx /. t -> 0
(*Aproximamos la segunda CI*) 
U[s_,1]:= U[s,1] =N[U[s,0]+Δt*Ut[s,0]+((((Δt))^2)/2)*Utt[s,0]];

(*Condiciones de contorno*)
U[0,r_]:=0; U[M+1,r_]:=0;
Ux[0,r_]:=0;Ux[M+1,r_]:=0;
Uxx[0,r_]:=0;Uxx[M+1,r_]:=0;

U[-1,r_]:=0; U[M+2,r_]:=0;
Ux[-1,r_]:=0;Ux[M+2,r_]:=0;
Uxx[-1,r_]:=0;Uxx[M+2,r_]:=0;

(*Definimos 3.3 y 3.4*)
Ux[s_,0]:=Ux[s,0]=(1/(12*Δt))(-U[s+2,0]+8*U[s+1,0]-8*U[s-1,0]+U[s-2,0]);
Uxx[s_,0]:=Uxx[s,0]=(1/(12*(((Δt))^2)))(-U[s+2,0]+16*U[s+1,0]-30*U[s,0]+16*U[s-1,0]-U[s-2,0]);
Ux[s_,1]:=Ux[s,1]=(1/(12*Δt))(-U[s+2,1]+8*U[s+1,1]-8*U[s-1,1]+U[s-2,1]);
Uxx[s_,1]:=Uxx[s,1]=(1/(12*(((Δt))^2)))(-U[s+2,1]+16*U[s+1,1]-30*U[s,1]+16*U[s-1,1]-U[s-2,1]);

(*Definimos los lambdas del sistema 3.6*)
l1[s_,r_]:=l1[s,r]=1-2*((α1)^2)*(((Δx))^2)*Uxx[s,r]-2*α1*α2*(((Δx))^2)*Uxx[s,r-1]-2*α1*α3*(((Δx))^2)*Uxx[s,r-2];
l2[s_,r_]:=l2[s,r]=-2*((α1)^2)*(((Δx))^2)*Ux[s,r]-4*α1*α2*(((Δx))^2)*Ux[s,r-1]-4*α1*α3*(((Δx))^2)*Ux[s,r-2];
l3[s_,r_]:=l3[s,r]=-1-α1*(((Δx))^2)-2*α1*α2*(((Δx))^2)*U[s,r-1]-2*α1*α3*(((Δx))^2)*U[s,r-2];
l4[s_,r_]:=l4[s,r]=(2+2*((α2)^2)*(((Δx))^2)*Uxx[s,r-1]+2*α2*α3*(((Δx))^2)*Uxx[s,r-2])*U[s,r-1]+
(-1+2*α2*α3*(((Δx))^2)*Uxx[s,r-1]+2*((α3)^2)*(((Δx))^2)*Uxx[s,r-2])*U[s,r-2]+
(2*((α2)^2)*(((Δx))^2)*Ux[s,r-1]+4*α2*α3*(((Δx))^2)*Ux[s,r-2])*Ux[s,r-1]+
2*((α3)^2)*(((Δx))^2)*Ux[s,r-2]*Ux[s,r-2]+(α2*(((Δx))^2)-2)*Uxx[s,r-1]+(α3*(((Δx))^2)+1)*Uxx[s,r-2];

W4[r_]:=W4[r]=Table[l4[i,r],{i,1,M}];

(*Definimos la matriz para resolver el sistema 3.6*)

T[r_]:=T[r]=SparseArray[{i_,j_}/;i==j+2->(((l2[i,r])/(12*(Δt)))-((l3[i,r])/(12*(((Δt))^2)))),{M,M}]+
SparseArray[{i_,j_}/;i==j+1->(-((2*l2[i,r])/(3*(Δt)))+((4*l3[i,r])/(3*(((Δt))^2)))),{M,M}]+
SparseArray[{i_,j_}/;i==j->(l1[i,r]-((5*l3[i,r])/(2*(((Δt))^2)))),{M,M}]+
SparseArray[{i_,j_}/;i==j-1->(((2*l2[i,r])/(3*(Δt)))+((4*l3[i,r])/(3*(((Δt))^2)))),{M,M}]+
SparseArray[{i_,j_}/;i==j-2->(-((l2[i,r])/(12*(Δt)))-((l3[i,r])/(12*(((Δt))^2)))),{M,M}];

(*Proceso de refinamiento iterativo para obtener (ux)r+1,s y (uxx)r+1,s*)
refinamiento[N0_,TOL_,M_, Ux0_,Ux1_,AuxUx2_, Uxx0_,Uxx1_,AuxUxx2_, α1_,α2_,α3_, iter_]:=
Module[{N0var= N0, Tolvar=TOL,Mvar=M,Ux0var = Ux0,Ux1var = Ux1,AuxUx2var = AuxUx2,Uxx0var = Uxx0,Uxx1var = Uxx1,AuxUxx2var = AuxUxx2,α1var=α1,α2var =α2,α3var=α3, itervar=iter,i,s,XX,X,Δx, Sol, UU, j,AuxUUx2var,AuxUUxx2var,AuxMat, AuxVect1, AuxVect2, AuxVect,l1,l2,l1new,l2new},
Δx = 200/(Mvar+1);
i=1;
XX=Table[X[s],{s,1, 2 M}];
l1=1-2*((α1)^2)*(((Δx))^2)*Uxx1var-2*α1*α2*(((Δx))^2)*Uxx1var-2*α1*α3*(((Δx))^2)*Uxx0var;
l2=-2*((α1)^2)*(((Δx))^2)*Ux1var-4*α1*α2*(((Δx))^2)*Ux1var-4*α1*α3*(((Δx))^2)*Ux0var;

AuxMat := AuxMat= SparseArray[{i_,j_}/;i==j->(-2*((α1var)^2)*(((Δx))^2)),{2 M,2 M}];
AuxVect1:= AuxVect1= 1-2*α1var*α2var*(((Δx))^2)*Uxx1var-2*α1var*α3var*(((Δx))^2)*Uxx0var-l1;

AuxVect2:= AuxVect2= -4*α1var*α2var*(((Δx))^2)*Ux1var-4*α1var*α3var*(((Δx))^2)*Ux0var - l2;

AuxVect=Join[AuxVect1, AuxVect2];

While[i≤N0var,
itervar=i;
Print["itervar = ",itervar];
Sol=NSolve[{AuxMat.XX==AuxVect},XX];
UU= XX/. Sol;

For[j=1,j≤ M,j++, AuxUUxx2var[j]= UU[[1,j]]];
For[j=1,j≤ M,j++, AuxUUx2var[j]= UU[[1,j+M]]];
AuxUx2var=Table[AuxUUx2var[s],{s,1,M}];
AuxUxx2var=Table[AuxUUxx2var[s],{s,1,M}];

If[Norm[AuxUxx2var-Uxx1var]<Tolvar && Norm[AuxUx2var-Ux1var]<Tolvar ,
{Print[" Control  YES"];
Return[AuxUx2var,Module];
Return[AuxUxx2var,Module];
Return[itervar ,Module]};
];

Print[" Control  NO"];

l1new=1-2*((α1)^2)*(((Δx))^2)*AuxUxx2var-2*α1*α2*(((Δx))^2)*Uxx1var-2*α1*α3*(((Δx))^2)*Uxx0var;
l2new=-2*((α1)^2)*(((Δx))^2)*AuxUx2var-4*α1*α2*(((Δx))^2)*Ux1var-4*α1*α3*(((Δx))^2)*Ux0var;


AuxVect1:= AuxVect1= 1-2*α1var*α2var*(((Δx))^2)*Uxx1var-2*α1var*α3var*(((Δx))^2)*Uxx0var-l1new;

AuxVect2:= AuxVect2= -4*α1var*α2var*(((Δx))^2)*Ux1var-4*α1var*α3var*(((Δx))^2)*Ux0var - l2new;

AuxVect=Join[AuxVect1, AuxVect2];

i++];

Return[AuxUx2var,Module];
Return[AuxUxx2var,Module];
Return[itervar ,Module];

]

N0=50;
TOL=0.0000001;
XX=Table[X[s],{s,1, M}];

(*Hacemos dos bucles para resolver nuestro sistema para cada valor de s y cada valor de r*)
For[r=2,r≤niteraciones,r++,     (*AQUÍ   EMPIEZA EL BUCLE TEMPORAL*)
{Print[" r = ", r],
iter=0;

For[s=1, s≤ M,s++,{
UU0[s]= U[s,r-2];
UUx0[s]= Ux[s,r-2];
UUxx0[s]= Uxx[s,r-2];

UU1[s]= U[s,r-1];
UUx1[s]= Ux[s,r-1];
UUxx1[s]= Uxx[s,r-1];

UU2[s]= U[s,r-1];
UUx2[s]= Ux[s,r-1];
UUxx2[s]= Uxx[s,r-1];
}];

UUx2[-1]:=0;UUx2[M+2]:=0;
UUxx2[-1]:=0;UUxx2[M+2]:=0;

UUUx0=Table[Ux[i,r-2],{i,1,M}];
UUUxx0=Table[Uxx[i,r-2],{i,1,M}];
UUUx1=Table[Ux[i,r-1],{i,1,M}];
UUUxx1=Table[Uxx[i,r-1],{i,1,M}];
UUUx2=Table[Ux[i,r-1],{i,1,M}];
UUUxx2=Table[Uxx[i,r-1],{i,1,M}];

refinamiento[N0,TOL,M,UUUx0,UUUx1,UUUx2, UUUxx0,UUUxx1,UUUxx2,α1,α2,α3, iter];


For[j=1,j≤ M,j++, Ux[j,r]=UUUx2[[j]]];
For[j=1,j≤ M,j++, Uxx[j,r]=UUUxx2[[j]]];


W4[r]=Table[l4[i,r],{i,1,M}];
UU:=UU=Table[0,{i,1,M}];
AA= T[r];
bb= W4[r];


XX=Table[X[s],{s,1, M}];
Sol=Solve[AA.XX== bb,XX];
UU= XX/. Sol;

For[s=1,s≤ M,s++, U[s,r]= UU[[1,s]]];

For[s=1,s≤ M,s++, Ux[s,r]=(1/(12*Δt))(-U[s+2,r]+8*U[s+1,r]-8*U[s-1,r]+U[s-2,r])];
For[s=1,s≤ M,s++, Uxx[s,r]=(1/(12*(((Δt))^2)))(-U[s+2,r]+16*U[s+1,r]-30*U[s,r]+16*U[s-1,r]-U[s-2,r])];

}
 ];    (*AQUÍ TERMINA EL BUCLE TEMPORAL*)

Plot1=Table[ListPlot[Table[{Δx s, U[s,r]},{s,0,M}], Joined-> True, PlotRange-> All,PlotStyle->Cyan, AxesLabel -> {"x", " "}, PlotLabel -> "T[x,t], t=" <> ToString[Δt r]], {r,0,niteraciones,25}]

Data3D = Table[{Δx s,0 , U[s,0]},{s,0,M+1}];
Do[Data3D=Join[Data3D,Table[{Δx s,Δt r , U[s,r]},{s,0,M+1}]],{r,1,niteraciones,10}]; 
Plot2=ListPlot3D[Data3D,PlotRange-> All,Mesh->{0,0},AxesLabel->{"Space - x","Time - t","u"},PlotLabel ->"Boussinesq improved equation",PlotStyle->Cyan]

t1=TimeUsed[ ];
Print[" El tiempo total empleado son ", t1-t0, " segundos."]