Timestep stability in a 1D heat diffusion model

Hi,

I have a 1D heat diffusion code which I was using on a timescale of 10s of years and I am now trying to use the same code to work on a scale of millions of years. Obviously if I keep my timestep the same this will take ages to calculate but if I increase my timestep I encounter numerical stability issues.

My questions are:

How should I approach this problem? What affects the maximum stable timestep? And how do I calculate this?

Many thanks,

Alex

close all
clear all

dx = 4;    % discretization step in m
dt = 0.0000001; % timestep in Myrs
h=1000;        % height of box in m
nx=h/dx+1;
model_lenth=1; %length of model in Myrs
nt=ceil(model_lenth/dt)+1;     % number of tsteps to reach end of model
kappa = 1e-6; % thermal diffusivity
x=0:dx:0+h;     % finite difference mesh
T=38+0.05.*x;  % initial T=Tm everywhere ...
time=zeros(1,nt);
t=0;
Tnew = zeros(1,nx);

%Lower sill
sill_1_thickness=18;
Sill_1_top_position=590;
Sill_1_top=ceil(Sill_1_top_position/dx);
Sill_1_bottom=ceil((Sill_1_top_position+sill_1_thickness)/dx);

%Upper sill
sill_2_thickness=10;
Sill_2_top_position=260;
Sill_2_top=ceil(Sill_2_top_position/dx);
Sill_2_bottom=ceil((Sill_2_top_position+sill_2_thickness)/dx);

%Temperature of dolerite intrusions
Tm=1300;

T(Sill_1_top:Sill_1_bottom)=Tm; %Apply temperature to intrusion 1

% unit conversion to SI:
secinmyr=24*3600*365*1000000;   % dt in sec
dt=dt*secinmyr;

%Plot initial conditions
figure(1), clf
f1 = figure(1); %Make full screen
set(f1,'Units', 'Normalized', 'OuterPosition', [0 0 1 1]); 
plot (T,x,'LineWidth',2)
xlabel('T [^oC]')
ylabel('x[m]')
axis([0 1310 0 1000])
title(' Initial Conditions')
set(gca,'YDir','reverse');

%Main calculation
for it=1:nt

      %Apply temperature to upper intrusion
      if it==10;  
          T(Sill_2_top:Sill_2_bottom)=Tm; 
      end

      for i = 2:nx-1
          Tnew(i) = T(i) + kappa*dt*(T(i+1) - 2*T(i) + T(i-1))/dx/dx;
      end

      Tnew(1) = T(1);
      Tnew(nx) = T(nx);

      time(it) = t;

      T = Tnew; %Set old Temp to = new temp for next loop
      tmyears=(t/secinmyr);

      %Plot a figure which updates in the loop of temperature against depth
      figure(2), clf
      plot (T,x,'LineWidth',2)
      xlabel('T [^oC]')
      ylabel('x[m]')
      title([' Temperature against Depth after ',num2str(tmyears),' Myrs'])
      axis([0 1300 0 1000])
      set(gca,'YDir','reverse');%Reverse y axis

      %Make full screen     
      f2 = figure(2); 
      set(f2,'Units', 'Normalized', 'OuterPosition', [0 0 1 1]); 

      drawnow

      t=t+dt;
  end