Definuj rychlostne vektorove pole

vp1[t_,x1_,x2_]=-(x2-0.5)
vp2[t_,x1_,x2_]=x1-0.5

Zaciatocna funkcia

phi0[x1_,x2_]=Exp[-100.0*((x1-.25)^2+(x2-0.5)^2)];
Plot3D[phi0[x,y],{x,0,1},{y,0,1},PlotRange->{0,1}]
ContourPlot[phi0[x,y],{x,0,1},{y,0,1},PlotRange->All]

Parametre pre numericku metodu, zaciatocna numericka funkcia

ii=40; L = 1.0; h=L/ii;
For[i=0, i<=ii, i++, 
For[j=0,j<=ii,j++,
phi[i,j,0] = phi0[i*h,j*h]]];
i=.;j=.;
ListPlot3D[Table[phi[i,j,0],{j,0,ii},{i,0,ii}],DataRange->{{0,1},{0,1}},PlotRange->All]
ListContourPlot[Table[phi[i,j,0],{j,0,ii},{i,0,ii}],DataRange->{{0,1},{0,1}},PlotRange->{0.0,1.0}]

Numericka metoda

T=Pi
tau=h/1.0
nn=IntegerPart[T/tau]
For[n=0,n<nn,n++,
{
f=Interpolation[Flatten[Table[{{i*h,j*h},phi[i,j,n]},{i,0,ii},{j,0,ii}],1],InterpolationOrder->2],
For[i=1, i<=ii-1, i++, 
For[j=1,j<=ii-1,j++,
{
v1=vp1[(n+1)*tau,i*h,j*h],
v2=vp2[(n+1)*tau,i*h,j*h],
x=i*h - tau/2*v1; y=j*h-tau/2*v2;
v1=vp1[(n+1/2)*tau,x,y],
v2=vp2[(n+1/2)*tau,x,y],
x=i*h- tau*v1; y=j*h-tau*v2;
phi[i,j,n+1]=f[x,y],
}
]
],
For[i=0, i<=ii, i++,phi[i,0,n+1]=phi[i,ii,n+1]=phi[0,i,n+1]=phi[ii,i,n+1]=0.0;],
If[ Mod[n,Quotient[nn,5]]==0,Print[ListPlot3D[Table[phi[i,j,n+1],{j,1,ii-1},{i,1,ii-1}],PlotRange->{0.0,1.0},DataRange->{{0,1},{0,1}}]]]
}
];
ListPlot3D[Table[phi[i,j,n],{j,1,ii-1},{i,1,ii-1}],PlotRange->{0.0,1.0}]

Animate[ListPlot3D[Table[phi[i,j,k],{j,1,ii-1},{i,1,ii-1}],PlotRange->{0.0,1.0}],{k,0,nn,Mod[nn,nn]+1},AnimationRepetitions->1]

Animate[ListContourPlot[Table[phi[i,j,k],{j,1,ii-1},{i,1,ii-1}]],{k,0,nn,1},AnimationRepetitions->1]

Min[Table[phi[i,j,nn],{j,1,ii-1},{i,1,ii-1}]]
Max[Table[phi[i,j,nn],{j,1,ii-1},{i,1,ii-1}]]