% This is an example of how to use command-line access
% to Matlab's PdeToolbox to solve a parabolic pde.
% The example is a variant of Matlab's, given in "help parabolic" or
% at http://www.mathworks.com/access/helpdesk/help/toolbox/pde/ug/parabolic.html.
% The comments are by DPO 04/2010.

% The problem:
%
%    The pde:
%         u_t - grad . ( grad u) = 2, inside the square -1 \le x,y \le 1.
%    The initial condition:
%         u(x,y,0) = 1 if (x,y) is in the circle 
%                         with radius 0.4 centered at the origin,
%         u(x,y,0) = 0 otherwise.                      
%    The boundary conditions:
%         u(x,y,t) = 0 if (x,y) is on the boundary of the square.
%    The times of interest:
%         0, 0.005, 0.010, ... , 0.100            

% Step 1: set up the x-y mesh.
%
% The mesh is specified by the 
%  --  Mesh points, stored in p.
%       (A mesh with n points makes p 2xn.  The x-coordinates
%        are in row 1 and the y-coordinates are in row 2.)
%  --  Triangle edge information, stored in e.
%       (This array is not so easy to decode.)
%  --  Triangle vertex information, stored in t.
%       (If column 10 of t is [212, 43, 5, 1]', then
%        the 10th triangle has mesh points 212, 43, and 5
%        as its vertices.)

% In this example, we let Matlab divide a square region into a
% coarse mesh of triangles, and then refine the mesh once.

[p,e,t] = initmesh('squareg');
[p,e,t] = refinemesh('squareg',p,e,t);

% Step 2: Specify initial conditions.
%
% Now we set the initial conditions v(x1,x2) for our pde. 
% In this example, we set v(x1,x2) = 1 if (x1,x2) is
% inside the circle with radius 0.4 centered at the origin,
% and v(x1,x2) = 0 otherwise.

u0 = zeros(size(p,2),1);
ix = find(sqrt(p(1,:).^2 + p(2,:).^2) < 0.4);
u0(ix) = 1;

% Step 3: Set the times of interest.
%
% Now we specify the times at which we want to compute
% the solution.

tlist = linspace(0, 0.1, 21);

% Step 4: Specify the PDE and the boundary conditions and solve.
%
% Now we tell Matlab to solve our problem.
% The boundary conditions are u(t,x) = 0 if x is on the
%     boundary of the square.  This is specified
%     by 'squareb1'.
%
% The last four arguments are  c, a, f, and d.  These
%     specify the pde to be
%      d u_t - grad . (c grad u) + a u = f.

u1 = parabolic(u0, tlist, 'squareb1', p, e, t, 1, 0, 2, 1);

% Step 5: Use the solution.
%
% The next command plots the solution at time 0.01.

pdeplot(p,e,t,'xydata',u1(:,3))
