% Parabolic Partial Differential Equations - bvp4c
% Parabolic partial differential equations (PDEs) are a class of PDE that is often used to model one-dimensional, time dependent problems. One example is one-dimensional, transient heat conduction. In this file I’ll describe how to solve simple parabolic PDEs using a built-in Matlab routine.
% To demonstrate the pdepe routine in Matlab, I’ll solve a 1-D transient conduction problem. The problem models a surface heat flux (q”) on an infinite slab with thermal conductivity k, density  \rho, and heat capacity c_p. The slab is of thickness L and the back surface is held at a temperature of 0 degrees for all time. The entire slab is initially at a temperature of 0 degrees and then the surface heat is turned on and held constant. The slab temperatures will rise until a steady state is reached. The equation we’re solving is
% \begin{array}{l}
% \rho c_p \frac{\partial T}{\partial t}=k\frac{\partial^2 T}{\partial x^2 }\\
% T\left(x,t=0\right)=0\\
% -k\frac{\partial T}{\partial x}=\textrm{q"}\;\textrm{at}\;x=0\\
% T=0\;\textrm{at}\;x=L
% \end{array}
% Matlab includes the routine pdepe to help solve parabolic PDEs. This routine requires the definition of several functions. As defined below, the pdex1pde function defines the partial differential equation itself. It is set up assuming that you are solving
% c\left(x,t,u,\frac{\partial u}{\partial x}\right)\frac{\partial u}{\partial t}=x^{-m} \frac{\partial }{\partial x}\left(x^m f\left(x,t,u,\frac{\partial u}{\partial x}\right)\right)+s\left(x,t,u,\frac{\partial u}{\partial x}\right)
% Here m=0 for Cartesian problems, 1 for cylindrical, and 2 for spherical. For our problem,
%  \begin{array}{l}
% c=\rho c_p \\
% f=k\frac{\partial T}{\partial x}\\
% s=0
% \end{array}
% The second function to be define is pdex1ic, which defines the initial conditions for the problem. In our case, we just define u_0 =0.
% The final function we require is pdex1bc, which defines the boundary conditions. The convention used in this routine is the boundary conditions are defined as follows:
% p\left(x,t,u\right)+q\left(x,t\right)\;f\left(x,t,u,\frac{\partial u}{\partial x}\right)=0
% We need to define this at both boundaries, so we define pl and ql on the left boundary and pr and qr on the right boundary. In our case, the right boundary does not involve a derivative of u, so we have qr=0 and pr=ur, implying that ur must be 0 for all time at the right boundary. On the left boundary, there is a derivative term, so we require pl to be set to the imposed heat flux and ql=1e6.
% The pdepe routine requires that we set up a series of mesh points for x and a series of time steps, and the linspace command is used for each of these. The plot command plots the temperature distribution after steady state has been reached.

global rho cp k
global q
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=10 %seconds
m = 0;
%Divide region 0<x<L into 1000 divisions
x = linspace(0,L,1000);
%Choose 100 time steps
t = linspace(0,tend,100);
sol = pdepe(m,@pdex1pde,@pdex1ic,@pdex1bc,x,t);
Temperature = sol(:,:,1);
plot(x,Temperature(end,:), 'LineWidth', 2)
xlabel('x (m)')
ylabel('temperature (C)')
title('Temperature Distribution at End of Simulation')
set(gca,'FontSize',12)
plot(t,Temperature(:,1), 'LineWidth', 2)
xlabel('time (s)')
ylabel('temperature (C)')
title('Temperature History of Surface')
set(gca,'FontSize',12)

%
% The PDE is set up here in pdex1pde.

function [c,f,s] = pdex1pde(x,t,u,DuDx)
global rho cp k
c = rho*cp;
f = k*DuDx;
s = 0;
end

%  The initial conditions are set up here in pdex1ic.

function u0 = pdex1ic(x)
u0 = 0;
end

%  And, finally, the boundaryt conditions are set up here in pdex1bc.

function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t)
global q
pl = q; 
ql = 1;
pr = ur;
qr = 0;
end

% © James Blanchard, University of Wisconsin-Madison, blanchard@engr.wisc.edu