function parabolic
global rho cp k
global q hcoef Tbulk
L=0.1     %m
k=200     %W/m-K
rho=10000 %kg/m^3
cp=500    %J/kg-K
q=1e6   %W/m^2
tend=100 %seconds
hcoef=1e4 %W/m^2-K
Tbulk=20  %Centigrade

m = 0;
x = linspace(0,L,200);
t = linspace(0,tend,50);

sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
% Extract the first solution component as u.
Temperature = sol(:,:,1);

% A solution profile can also be illuminating.
figure, plot(x,Temperature(end,:))
title(strcat('Solution at t = ', num2str(tend)))
xlabel('Distance x')
ylabel('Temperature (C)')

%Plot surface temperature vs. time
figure, plot(t,Temperature(:,1))
title('Surface Temperature')
xlabel('Time (s)')
ylabel('Temperature (C)')
% --------------------------------------------------------------
function [c,f,s] = pdex1pde(x,t,u,DuDx)
global rho cp k
c = rho*cp;
f = k*DuDx;
s = 0;
% --------------------------------------------------------------
function u0 = pdex1ic(x)
global q hcoef Tbulk
u0 = Tbulk;
% --------------------------------------------------------------
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
global q hcoef Tbulk
pl = q;
ql = 1;
pr = hcoef*(ur-Tbulk);
qr = 1;
