function pdex3
%PDEX3  Example 3 for PDEPE
%   This example illustrates the use of PDEVAL to obtain partial derivatives
%   that are part of solving the problem. A relatively fine spatial mesh is
%   needed to obtain accurate partial derivatives. Values are required at a
%   relatively large number of times.
%
%   This problem arises in transistor theory. It is used in [1] to illustrate
%   solution of PDEs with series. In the form expected by PDEPE, the single PDE
%   is
%
%    [1] .*  D_ [u] =  D_ [   d*Du/Dx   ] +  [ -(eta/L)*Du/Dx ]
%            Dt        Dx
%    ---        ---        -------------      ----------------
%     c          u         f(x,t,u,Du/Dx)      s(x,t,u,Du/Dx)
%
%   Here d and eta are physical constants. The equation is to hold on an
%   interval 0 <= x <= L for times t >= 0. For the problem at hand, L = 1.
%   The initial condition
%
%      u(x,0) = (K*L/d)*((1 - exp(-eta*(1 - x/L)))/eta)
%
%   involves another physical constant K. The left bc is u(0,t) = 0:
%
%      [u]    +     [0] .* [       Du/Dx          ] = [0]
%
%      ---          ---    ------------------------   ---
%    p(0,t,u)      q(0,t)        f(0,t,u,Du/Dx)        0
%
%   The right bc is u(L,t) = 0:
%
%      [u]    +     [0] .* [       Du/Dx          ] = [0]
%
%      ---          ---    ------------------------   ---
%    p(L,t,u)      q(L,t)        f(L,t,u,Du/Dx)        0
%
%   See the nested functions PDEX3PDE, PDEX3IC, and PDEX3BC for the coding of
%   the problem definition.
%
%   A quantity of physical interest is the emitter discharge current
%
%      I(t) = (I_p*d/K)*Du(0,t)/Dx
%
%   where I_p is another physical constant. This is meaningful only for t > 0
%   because of the inconsistency in the boundary values at x = 0 for t = 0 and
%   t > 0. This is reflected in the failure to converge of the series solution
%   of [1] for t = 0. Because values for the physical parameters are not
%   specified in [1], nominal values are assumed. I(t) changes rapidly, so
%   output at a relatively large number of t are needed. Recall that the number
%   and placement of the entries of t have little effect on the cost of the
%   computation. Also, the approximation of the solution is only second order
%   in space and the approximation of the partial derivative is of still lower
%   order. Nevertheless, there is reasonable agreement between the computed
%   emitter discharge current and that obtained from the series.
%
%   [1] E.C. Zachmanoglou and D.L. Thoe,Introduction to Partial Differential
%       Equations with Applications, Dover, New York, 1986.
%
%   See also PDEPE, PDEVAL, FUNCTION_HANDLE.

%   Lawrence F. Shampine and Jacek Kierzenka
%   Copyright 1984-2014 The MathWorks, Inc.

% Problem parameters, shared with nested functions.
L = 1;
d = 0.1;
eta = 10;
K = d;
I_p = 1;

m = 0;
x = linspace(0,L,41);
t = linspace(0,1,51);
sol = pdepe(m,@pdex3pde,@pdex3ic,@pdex3bc,x,t);
u = sol(:,:,1);

figure;
surf(x,t,u);
title('Numerical solution computed with 41 mesh points.');
xlabel('Distance x');
ylabel('Time t');

nt = length(t);
I = zeros(1,nt);
seriesI = zeros(1,nt);
iok = 2:nt;
for j = iok
   % At time t(j), compute Du/Dx at x = 0.
   [~,I(j)] = pdeval(m,x,u(j,:),0);
   seriesI(j) = serex3(t(j));
end
% I(t) = (I_p*d/K)*Du(0,t)/Dx
I = (I_p*d/K)*I;

figure;
plot(t(iok),I(iok),t(iok),seriesI(iok));
legend('From PDEPE','From series');
title('Emitter discharge current I(t)');
xlabel('Time t');

% -----------------------------------------------------------------------
% Nested functions -- parameters are provided by the outer function.
%

   function [c,f,s] = pdex3pde(x,t,u,DuDx)
      % Evaluate the differential equations components.
      c = 1;
      f = d*DuDx;
      s = -d*eta*DuDx;
   end
% -----------------------------------------------------------------------

   function u0 = pdex3ic(x)
      % Evaluate the initial conditions components.
      u0 = (K*L/d)*(1 - exp(-eta*(1 - x)))/eta;
   end
% -----------------------------------------------------------------------

   function [pl,ql,pr,qr] = pdex3bc(xl,ul,xr,ur,t)
      % Evaluate the boundary conditions components.
      pl = ul;
      ql = 0;
      pr = ur;
      qr = 0;
   end
% -----------------------------------------------------------------------

   function It = serex3(t)
      % Approximate I(t) by 40 terms of a series expansion.
      It = 0;
      for n = 1:40
         temp = (n*pi)^2 + 0.25*eta^2;
         It = It + ((n*pi)^2 / temp)* exp(-(d/L^2)*temp*t);
      end
      It = 2*I_p*((1 - exp(-eta))/eta)*It;
   end
% -----------------------------------------------------------------------

end  % pdex3
