clear
clf;
nit = 10000;
f = 1;
d = 5;
N = 50;
dx = N/5*d;
dt = 1000000/24;
cp = 1;
k = .03;
rho = 2.7;
c = k/(cp*rho);
eta = dt*c/ dx^2;
for i = [2:(1/5)*N]
    T(i) = 0;
end
for i = [(1/5)*N:(2/5)*N]
    T(i) = 1000;
end
for i = [(2/5)*N:(3/5)*N]
    T(i) = 0;
end
for i = [(3/5)*N:(4/5)*N]
    T(i) = 800;
end
for i = [(4/5)*N:N+1]
    T(i) = 0;
end
plot(T)
for ii = 1:nit
    ii;
    Tnew = T;
    for i = 2:N
        Tnew(i) = T(i) + eta* (T(i+1) - 2.0* T(i) + T(i-1));
    end
    T = Tnew;
    if mod(ii,100)==0 
       plot(T,'-ro')
       axis([0 50 0 1000])
       drawnow
       ii;
       saveas(gcf, ['lava',sprintf('%03d',f),'.png'])
       f = f+1;
    end

end