function [ x, t, u ] = method_of_lines ( xrange, trange, nx, nt )
%
%  function [ x, t, u ] = method_of_lines ( xrange, trange, nx, nt )
%
%  Purpose:
%
%    METHOD_OF_LINES solves a simple PDE using the method of lines.
%
%  Discussion:
%
%    The PDE has the form:
%
%      Ut = Uxx + G(X,T) 
%
%    for
% 
%      XMIN < X < XMIN, 
%      0 <= T
%
%    with an initial condition of:
%
%      U(X,0) = SIN ( PI * X )
%
%    and boundary conditions
%
%      U(XMIN,T) = 0
%      U(XMAX,T) = 0
%
%    The algorithm uses a simple Euler method in time, after discretizing 
%    the second order space derivative using central differences.
%
%    The function G(X,T) is "cooked up" to be
%
%      G(X,T) = ( pi*pi - 1/10 ) * exp ( - T / 10 ) * sin ( pi * X )
%
%    so that the exact solution of the PDE is
%
%      U(X,T) = exp ( - T / 10 ) * sin ( pi * X )
%
%  Modified:
%
%    04 February 2000
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    XRANGE(2) is the minimum and maximum X value;
%    TRANGE(2) is the minimum and maximum T value;
%    NX is the number of X intervals (NX+1 points);
%    NT is the number of T intervals (NT+1 time steps).
%
%  Output:
%
%    X(NX+1), the X points;
%    T(NT+1), the T points;
%    U(NX+1,NT+1), the computed solution.
%

%
%  Initialize.
%
%  The spatial domain is [XMIN,XMAX], so NX defines the value of DX.
%
xmin = xrange(1);
xmax = xrange(2);
dx = ( xmax - xmin ) / nx;
x = linspace ( xmin, xmax, nx+1 );
%
%  The time step is fixed at 0.002, so NT defines how far we advance in time.
%
tmin = trange(1);
tmax = trange(2);
dt = ( tmax - tmin ) / ( nt - 1 );
t = linspace ( tmin, tmax, nt+1 );
%
%  Define U to be an array of the right size, and all zero.
%
u = zeros ( nx+1, nt+1 );
%
%  Set the initial condition.
%
it = 1;
for ix = 1 : nx+1
  u(ix,it) = sin ( pi * x(ix) );
end
%
%  Take a time step.
%
%  The simple boundary conditions are handled by fixing the first and
%  last entries of U.  
%
%  The other values advance in time using the Euler method.
%  The second derivative Uxx is approximated by a central difference.
%  If we know the exact solution, then the function G(X,T) 
%  is simply (Ut-Uxx).
%
for it = 2 : nt+1

  u(1,it) = u(1,it-1);

  for ix = 2 : nx
    u(ix,it) = u(ix,it-1) + dt * (...
      ( u(ix+1,it-1) - 2 * u(ix,it-1) + u(ix-1,it-1) ) / dx^2 ...
      + ( pi*pi - 0.1 ) * exp ( - t(it-1) / 10.0 ) * sin ( pi * x(ix) ) );
  end

  u(nx+1,it) = u(nx+1,it-1);

end