
Get["equationDiffusion.m"];
n = 200;
u = Table[0,{n}];
s = Table[0,{n}];
coef = {{1,1}};
t = 0;
{u1,t}=diffusion[n,"dirichlet",1,"dirichlet",0,coef,s,u,t,0.0001,0.001];
{u2,t}=diffusion[n,"dirichlet",1,"dirichlet",0,coef,s,u1,t,0.001,0.01];
{u3,t}=diffusion[n,"dirichlet",1,"dirichlet",0,coef,s,u2,t,0.01,0.1];
{u4,t}=diffusion[n,"dirichlet",1,"dirichlet",0,coef,s,u3,t,0.1,1];
                
ListLinePlot[{u1,u2,u3,u4},PlotRange->{{0,n},{0,1}},AxesLabel->{"x","T"}]

n = 500;
u = Table[0,{n}];
s = Table[0,{n}];
For[j=1,j<=Floor[n/2],j++,u[[j]]=1];
coef = {{0.45,1},{0.55,0.05},{1,1}};
t = 0;
{u1,t}=diffusion[n,"neumann",0,"neumann",0,coef,s,u,t,0.00001,0.001];
{u2,t}=diffusion[n,"neumann",0,"neumann",0,coef,s,u1,t,0.001,0.01];
{u3,t}=diffusion[n,"neumann",0,"neumann",0,coef,s,u2,t,0.01,0.1];
{u4,t}=diffusion[n,"neumann",0,"neumann",0,coef,s,u3,t,0.1,1];
                
ListLinePlot[{u1,u2,u3,u4},PlotRange->{{0,n},{0,1}},AxesLabel->{"x","T"}]
