clf
clear
 
N =70;
L= 1500;

 
p=2.65; %density
k=(0.031)/(1.1*2.65);
x=0:L/(8.5*N):L;
y = 1:L; %
c=1.1;
 
T=zeros(size(x));
 
dx=L/(8*N);
dt= 310; 
T(1:N)=0;   %temp of host rock
T(N+1:1.5*N)=1300;  %temp of dike 
T(1.5*N+1:2.5*N)=0;  
T(2.5*N+1:3*N)=1200;
T(3*N+1:4*N)=0; 
T(4*N+1:5*N)=1100;
T(5*N+1:6*N)=0;
T(6*N+1:6.5*N)=1000;
T(6.5*N+1:7*N)=0;
T(7*N+1:7.5*N)=900;
T(7.5*N+1:8.52*N)=200;

f=0;
for t=1:3600 %time step in seconds
  Tnew=T;
  for i=2:8.5*N
      Tnew(i)=(dt*k/dx^2)*(T(i+1)-2*T(i)+T(i-1))+T(i);
  end
  T=Tnew;
  plot(x,T,'o-')
  grid on;
  axis([0 1500 0 1300]);
  pause(.0001);
  title('Lava Dikes cooling over time');
  ylabel('Temperature (C)'); 
  xlabel('Distance from origin (decimeters)');

 if mod(t,120)==0;
      f = f+1;
     %surf(x,y,T);
  saveas(gcf, ['final_dikes',sprintf('%03d',f),'.png']);
  if f>99;
      break
  end
 end
en