clear;
clf;
lambda=10;
rhoc=1380000;
a=lambda/rhoc;
q=0;
deltat=100;
tmax=86400;
initemp=283.15;
deltax=0.1;
jmax=10/deltax;
bc=50;
ts=0;
tempe=zeros(jmax,1);
for j=1:jmax
    temp(j)=initemp;
end
//    start calculation
for t=1:deltat:tmax //    time loop
    temp(1) = initemp+bc;       // setting Dirichlet bc
    tempe(jmax)=tempe(jmax-1);  // setting Neumann bc
    for i=2:(jmax-1) 
        //     now the main calculation
        tempe(i)=temp(i)+(a*deltat/deltax^2)*(temp(i-1)-2.0*temp(i)+temp(i+1))+q;        
    end //     end of grid points loop
    tempe(1) = temp(1);     // restoring Dirichlet bc
    tempe(jmax)=tempe(jmax-1); // setting Neumann bc
    //     delivery of the provisional result back to the time loop value
    for i=1:jmax
        temp(i)=tempe(i);
    end
    if(modulo((t-1),3600)==0)
         ts = ts + 1
         t
         for i=1:jmax
             tstore(i,ts)=temp(i);
         end
         plot(tstore);
     end
end     //     end of time loop
