pb=10;

A=SparseArray[
{
{i_,i_}->2,
{i_,j_}/;Abs[i-j]==1->-1
},{pb-1,pb-1}];

Identita=SparseArray[{i_,i_}->1,{pb-1,pb-1}];

L=1.0 ;
h=L/pb

grid=Range[h,L-h,h];

phi0[x_]:=Sin[1.0*Pi*x](*presna pociatocna podmienka*)

Plot[phi0[x],{x,0,L}]

phi[0]=phi0[grid](*0 je index pre cas   dava to vektor v priestore !!!*)

ListLinePlot[Transpose[{Prepend[Append[grid,L],0],Prepend[Append[phi[0],0.0],0.0]}]](*numericka nulta pociatocna podmienka*)

tend=0.1;pc=10;dt=tend/pc
alpha=0.5;
G=1.;

solver=LinearSolve[h/dt Identita+alpha G/h A];

For[
n=0,n<pc,n++,
{
b=h/dt Identita.phi[n]-(1.0-alpha)*G/h A.phi[n],
phi[n+1]=solver[b]
}
]

ListLinePlot[Transpose[{grid,phi[pc]}]]

Animate[ListLinePlot[Transpose[{grid,phi[n]}],PlotRange->{-1,1}],{n,0,pc,1},AnimationRunning->False]

(*pridame konvekciu !!! *)

U=SparseArray[
{
{i_,i_}->1,
{i_,j_}/;i-j==1->-1
},{pb-1,pb-1}];

MatrixForm[U]

V=SparseArray[
{
{i_,i_}->-1,
{i_,j_}/;i-j==-1->1
},{pb-1,pb-1}];

MatrixForm[V]

konvekcia[ pb_,L_]:=Module[{A,U,V},
A=SparseArray[
{
{i_,i_}->2,
{i_,j_}/;Abs[i-j]==1->-1
},{pb-1,pb-1}];
Print[MatrixForm[A]];
Identita=SparseArray[{i_,i_}->1,{pb-1,pb-1}];
U=SparseArray[
{
{i_,i_}->1,
{i_,j_}/;i-j==1->-1
},{pb-1,pb-1}];
V=SparseArray[
{
{i_,i_}->-1,
{i_,j_}/;i-j==-1->1
},{pb-1,pb-1}];
Print[MatrixForm[U]];
Print[MatrixForm[V]];


h=L/pb;
grid=Range[h,L-h,h];
phi0[x_]:=Sin[1.0*Pi*x];
Plot[phi0[x],{x,0,L}];
phi[0]=phi0[grid];


tend=0.1;pc=10;dt=tend/pc;
alpha=0.5;
beta=0.5;
G=1.;

solver=LinearSolve[h/dt Identita+alpha G/h A];

For[
n=0,n<pc,n++,
{
b=h/dt Identita.phi[n]-(1.0-alpha)*G/h A.phi[n],
phi[n+1]=solver[b]
}
];

a=ListLinePlot[Transpose[{grid,phi[pc]}]];
aa=Animate[ListLinePlot[Transpose[{grid,phi[n]}],PlotRange->{-1,1}],{n,0,pc,1},AnimationRunning->False]

]

konvekcia[10,1.0]

ll

pb=100;L=2.5;

ρ=1;
u=-100;

A=SparseArray[
{
{i_,i_}->2,
{i_,j_}/;Abs[i-j]==1->-1
},{pb-1,pb-1}];

Identita=SparseArray[{i_,i_}->1,{pb-1,pb-1}];
U=SparseArray[
{
{i_,i_}->1,
{i_,j_}/;i-j==1->-1
},{pb-1,pb-1}];
V=SparseArray[
{
{i_,i_}->-1,
{i_,j_}/;i-j==-1->1
},{pb-1,pb-1}];
(*Print[MatrixForm[A]];
Print[MatrixForm[U]];
Print[MatrixForm[V]];
*)

h=L/pb;
grid=Range[h,L-h,h];
phi0[x_]:=Sin[1.0*Pi*x];

σ=100;
xc=5;
(*phi0[x_]:=E^(-6*(x-xc)^2);*)


phi[0]=phi0[grid];


tend=2;pc=100;dt=tend/pc;
alpha=0.5;
beta=0.5;
G=1.;

solver=LinearSolve[h/dt Identita+alpha G/h A+beta*(Max[ρ*u,0]U+Min[ρ*u,0]V)];

For[
n=0,n<pc,n++,
{
b=h/dt Identita.phi[n]-(1.0-alpha)*G/h A.phi[n]+(1-beta)*(Max[ρ*u,0]U.phi[n]+Min[ρ*u,0]V.phi[n]),
phi[n+1]=solver[b]
}
];

ListLinePlot[Transpose[{grid,phi[pc]}]]
riesenie[x_,t_]:=(1/Sqrt[1+4*σ*t])*E^(-σ*((x-xc-ρ*u*t)^2/(1+4*σ*t)));
(*Animate[{Plot[riesenie[x,n],{x,0,L}],ListLinePlot[Transpose[{grid,phi[n]}],PlotRange->{-1,1}]},{n,0,pc,1},AnimationRunning->False]*)
Animate[ListLinePlot[Transpose[{grid,phi[n]}],PlotRange->{-1,1}],{n,0,pc,1},AnimationRunning->False]

.

.

σ=100;
xc=5;
t=1;
riesenie[x_,t_]:=(1/Sqrt[1+4*σ*t])*E^(-σ*((x-xc-ρ*u*t)^2/(1+4*σ*t)));
Animate[Plot[riesenie[x,t],{x,0,20}],{t,0,5},AnimationRunning->False]

pb=100;L=1.0;

ρ=1;
u=1;

phiL=10;


A=SparseArray[
{
{i_,i_}->2,
{i_,j_}/;Abs[i-j]==1->-1
},{pb-1,pb-1}];

Identita=SparseArray[{i_,i_}->1,{pb-1,pb-1}];
U=SparseArray[
{
{i_,i_}->1,
{i_,j_}/;i-j==1->-1
},{pb-1,pb-1}];
V=SparseArray[
{
{i_,i_}->-1,
{i_,j_}/;i-j==-1->1
},{pb-1,pb-1}];
(*Print[MatrixForm[A]];
Print[MatrixForm[U]];
Print[MatrixForm[V]];
*)

h=L/pb;
grid=Range[h,L-h,h];
(*phi0[x_]:=Sin[1.0*Pi*x];*)

σ=100;
xc=5;
(*phi0[x_]:=E^(-6*(x-xc)^2);*)
phi0[x_]:=Sin[1.0*Pi*x];

phi[0]=phi0[grid];


tend=2;pc=100;dt=tend/pc;
alpha=0.5;
beta=0.5;
G=1.;

solver=LinearSolve[h/dt Identita+alpha G/h A+beta*(Max[ρ*u,0]U+Min[ρ*u,0]V)];

For[
n=0,n<pc,n++,
{
b=h/dt Identita.phi[n]-(1.0-alpha)*G/h A.phi[n]+(1-beta)*(Max[ρ*u,0]U.phi[n]+Min[ρ*u,0]V.phi[n]),
phi[n+1]=solver[b]
}
];

ListLinePlot[Transpose[{grid,phi[pc]}]]
riesenie[x_,t_]:=(1/Sqrt[1+4*σ*t])*E^(-σ*((x-xc-ρ*u*t)^2/(1+4*σ*t)));
(*Animate[{Plot[riesenie[x,n],{x,0,L}],ListLinePlot[Transpose[{grid,phi[n]}],PlotRange->{-1,1}]},{n,0,pc,1},AnimationRunning->False]*)
Animate[ListLinePlot[Transpose[{grid,phi[n]}],PlotRange->{-1,1}],{n,0,pc,1},AnimationRunning->False]