(* parabolicke rovnice *)
f[x_,t_]:=2*t*x-1;
g[x_]:=x^2;
beta[t_]:=1+t+t^2;
alfa[t_]:=t;
pres[x_,t_]:=x*t^2+x^2+t;
a=0;
b=1;
k=1;
pn=6;
cn=5;

(* vypocet pociatocneho riesenia *)
h=N[(b-a)/(pn+1)]
tau=h^2/(2*k)

(* vypocet, grafy *)

xs=Array[1,pn+1,0];
up=Array[1,pn+1,0];
u=Array[1,pn+1,0];
Do[
{xs[[i]]=a+i*h,up[[i]]=g[a+i*h]},{i,0,pn+1}
];
numer[0]=Table[{0,h*i1,up[[i1]]},{i1,0,pn+1}];
num=Table[{xs[[i]],up[[i]]},{i,0,pn+1}];
numgraf=ListPlot[num,PlotStyle->{PointSize[0.02],RGBColor[1,0,1]}];
presgraf=Plot[pres[x,0],{x,a,b}];
Show[presgraf,numgraf]

con=(tau*k)/(h^2);
tlac=1;
Print[" kontrola stability: ",con];

Do[
{
Do[
{
u[[i]]=up[[i]]+con*(up[[i-1]]-2*up[[i]]+up[[i+1]])+tau*f[a+i*h,(j-1)*tau]
},{i,1,pn}
]
Do[up[[i]]=u[[i]],{i,1,pn}
],
up[[0]]=alfa[j*tau],
up[[pn+1]]=beta[j*tau],
num=Table[{xs[[i1]],up[[i1]]},{i1,0,pn+1}],
numer[j]=Table[{tau*j,h*i1,up[[i1]]},{i1,0,pn+1}],
If[Mod[j,tlac]==0,{numgraf=ListPlot[num,PlotRange->{{0,1},{0,1.5}},PlotStyle->{PointSize[0.02],RGBColor[1,0,1]}],
presgraf=Plot[pres[x,j*tau],{x,a,b},PlotRange->{{0,1},{0,1.5}}],
Print["cas ", j*tau],
spolu=Show[presgraf,numgraf],
Print[spolu]
}]
},{j,1,cn}
]

a=Table[numer[j],{j,0,cn}];
aha=Flatten[a,1];
aprp1=ListPlot3D[aha,Mesh->All]