heat.m
function heat_neil
% Numerical Modeling in the Geosciences, Humphrey
% homework problem in heat flow
% Mainly to get you working with Matrices
% ***************************************************************************
% Problem: steady state temperature field in a 1000 m thick layer of rock,
% with a surface temperature of 0degrees and a heat flux of 65 mW per square meter
% ***************************************************************************
clear
K = 2.2; % granite
n = 9; %number of internal nodes, total nodes (including BCs) =n+1
thickness = 1000; %thickness of rock layer in vertical, in meters
delx = thickness/(n+1);
% delx used for plotting
upperBC = 0; % temp BCs, in degrees
upperq = 6.5e-2; %
T(1) = upperBC;
T(2) = upperq*delx/K;
for j=3:n+2
T(j) = T(j-1) + upperq*delx/K;
end
%plot result (I always like to plot the results to see if there are mistakes
plot(T, 0:delx:thickness,'-o');
axis ij; % makes the upper left corner the origin
axis([0 30 0 1000]);
title('Steady State temperature in a rock slab');
ylabel('Depth from surface in meters');
xlabel('Temperature in degrees');