function [t,u]=pdefdm(npde,jmax,tspan,u0,ut0,PDE_F,PDE_G,BC_L,BC_R,a,b,reltol,abstol)
%
% [t,u]=pdefdm(npde,jmax,tspan,u0,ut0,PDE_F,PDE_G,BC_L,BC_R,a,b,reltol,abstol)
%
% This function integrates a system of (npde) one-dimensional PDEs in the form
%
%   F(t, x, u, u_x, u_t) = (d/dx) G (t, x, u, u_x, u_t)  for a < x < b
%                      0 = B_L(t, x, u, u_x, u_t)         at x = a
%                      0 = B_R(t, x, u, u_x, u_t)         at x = b
%
% using a fixed mesh.
%  
% The physical and mesh PDEs are discretized using central finite differences and
% the resulting ODE system is integrated using the Matlab built-in solver ode15i.
%
% Note that the first 11 arguments must be specified by the user. The other 2
% arguments (reltol,abstol) are optional.
%
% variables:
%
% npde:     the number of physical PDEs contained in the system to be integrated.
% jmax:     the number of mesh points.
% tspan:    an array specifying the time interval for integration.
% u0:       array of size npde*jmax, intial values for the physical variables
%           and the mesh.
%           u0(npde*(j-1)+l): j = 1:jmax, l=1:npde
%                             the initial value for the l-th component of the
%                             physical variable at the mesh point x_j.
% ut0:      array of size npde*jmax, approximate values for the initial time
%           derivatives of the physical variables and mesh. approximate
%           values for the time derivatives are required by the implicit DAE solver
%           ode15i.
% PDE_F,PDE_G,BC_L,BC_R: functions used to define the system of physical PDEs:
%
%           PDE_F := F(t,x,u,u_x,u_t)
%           PDE_G := G(t,x,u,u_x,u_t)
%            BC_L := B_L(t,x,u,u_x,u_t)
%            BC_R := B_R(t,x,u,u_x,u_t)
%
%    function f=PDE_F(t,x,u,ux,ut)
%    function g=PDE_G(t,x,u,ux,ut)
%
%           input variables:
%           t is scalar, x=x(1:m) is a row vector, and u=u(1:npde,1:m), ux=ux(1:npde,1:m),
%           and ut=ut(1:npde,1:m) are matrices.
%
%           output variable:
%           f=f(1:npde,1:m) and g=g(1:npde,1:m) should be matrices.
%
%    function res=BC_L(t,x,u,ux,ut)
%    function res=BC_R(t,x,u,ux,ut)
%
%           input variables:
%           t and x are scalar and u=u(1:npde), ux=ux(1:npde), and ut=ut(1:npde) are
%           column vectors.
%
%           output variable:
%           res=res(1:npde) should be a column vector.
%
% a, b:     the left and right ends of the physical interval.
% reltol, abstol: relative and absolute tolerances used in integration of
%           the ode system.
%

%
% Copyright (C) 2010 Weizhang Huang and Robert D. Russell
% all rights reserved.
%
% This program is provided "as is", without warranty of any kind.
% Permission is hereby granted, free of charge, to use this program
% for personal, research, and education purposes. Distribution or use
% of this program for any commercial purpose is permissible
% only by direct arrangement with the copyright owner.
%

% the other functions contained in this file:
%
% 1. res=pdefdm_res(t,v,vt,PDE_F,PDE_G,BC_L,BC_R,jmax,npde,TauFun,monitor)
% 2. [SPDY,SPDYP] = jpattern(npde,jmax)
%

   % set the default values for the parameters

   if nargin < 11
       error('The first 11 arguments must be defined.');
   elseif nargin < 12
      reltol=[];
      abstol=[];
   elseif nargin < 13
      abstol=[];
   elseif nargin > 13
      error('Two many arguments are specified.');
   end    

   flag_reltol = 0;
   flag_abstol = 0;
   if isempty(reltol) | (reltol<=0)
      reltol=1e-6;
      flag_reltol = 1;
   end
   if isempty(abstol) | (abstol<=0)
      abstol=1e-4;
      flag_abstol=1;
   end

% set options

   nn = npde*jmax;
			   
   [SPDY,SPDYP]=jpattern(npde,jmax);
   j=round(jmax/2);
   if j < 1, j = 1; end
   j = npde*(j-1)+1;
   opts=odeset('RelTol',reltol,'AbsTol',abstol,'OutputFcn',@odeprint,...
               'OutputSel',j,'JPattern',{SPDY,SPDYP},'Stats','on');
      
% use a helper to find better initial derivatives

   free_y=zeros(nn,1);
   
   [u0md,ut0]=decic(@pdefdm_res,tspan(1),u0,free_y,ut0,[],opts,...
               PDE_F,PDE_G,BC_L,BC_R,a,b,jmax,npde);
   dd1 = norm(u0md(1:nn)-u0(1:nn),inf);
   dd2 = norm(u0md((nn+1):end)-u0((nn+1):end),inf);
   
% call ode15i

   [t,u]=ode15i(@pdefdm_res,tspan,u0md,ut0,opts,...
                PDE_F,PDE_G,BC_L,BC_R,a,b,jmax,npde);
                
% print messages

if flag_reltol==1
   fprintf('\n *** The default value reltol=1e-6 is used.\n');
end
if flag_abstol==1
   fprintf('\n *** The default value abstol=1e-4 is used.\n');
end
if (dd1+dd2)>eps
   fprintf('\n *** Warning: a modified initial solution is used.\n');
     fprintf('     max diff in solution = %e\n', dd1);
     fprintf('         max diff in mesh = %e\n', dd2);
end

% end of pdefdm

% --------------------------------------------------------------------------

function res=pdefdm_res(t,v,vt,PDE_F,PDE_G,BC_L,BC_R,a,b,jmax,npde)
%
% this function computes the res function to be called by an ode
% solver. the res function results from the central finite difference
% discretizetion of the (npde) physical pdes and the corresponding
% boundary conditions.
%

   % declare arrays
   
   Res=zeros(npde,jmax);
   
   u=zeros(npde,1);
   ux=zeros(npde,1);
   ut=zeros(npde,1);
   XX = zeros(1,jmax-1);
   U = zeros(npde,jmax-1);
   UX = zeros(npde,jmax-1);
   UT = zeros(npde,jmax-1);
   G = zeros(npde,jmax-1);
   nn = npde*jmax;
   
   % extract physical variables and mesh from v and vt

   V = reshape(v(1:nn),npde,jmax);
   Vt = reshape(vt(1:nn),npde,jmax);
   x = linspace(a,b,jmax);

   % for physical pdes: interior points
   
   % for F term
   for j=2:(jmax-1)
	  XX(j) = x(j);
      U(:,j) = V(:,j);
      UX(:,j) = (V(:,j+1)-V(:,j-1))/(x(j+1)-x(j-1));
      UT(:,j) = Vt(:,j);
   end
   % Res(:,1) will be reset by boundary conditions
   XX(1) = x(1);
   U(:,1) = V(:,1);
   UX(:,1) = (V(:,2)-V(:,1))/(x(2)-x(1));
   UT(:,1) = Vt(:,1);
   Res(:,1:(jmax-1)) = feval(PDE_F,t,XX,U,UX,UT);  
   
   % for G term
   for j=1:(jmax-1)
	  XX(j) = (x(j+1)+x(j))*0.5;
      U(:,j) = (V(:,j+1)+V(:,j))*0.5;
      UX(:,j) = (V(:,j+1)-V(:,j))/(x(j+1)-x(j));
	  UT(:,j) = 0.5*(Vt(:,j+1)+Vt(:,j));
   end
   G = feval(PDE_G,t,XX,U,UX,UT);
   for j=2:(jmax-1)
	  Res(:,j) = Res(:,j) - (G(:,j)-G(:,j-1))*2/(x(j+1)-x(j-1));
   end
   
   % for physical pdes: boundary points

   % at the left end
   j = 1;   
   xx = x(j);
   dx = 1.0/(x(j+1)-x(j));
   dx1 = x(2)-x(1);
   dx2 = x(3)-x(1);
   u = V(:,j);
   % 1st order approximation
%   ux = dx*(V(:,j+1)-V(:,j));
   % 2nd order approximation
   ux = (dx1*dx1*V(:,3)-dx2*dx2*V(:,2)-(dx1*dx1-dx2*dx2)*V(:,1)) ...
        /((dx1-dx2)*dx1*dx2);
   ut = Vt(:,j);
   Res(:,j) = feval(BC_L,t,xx,u,ux,ut);
   
   % at the right end
   j = jmax;   
   xx = x(j);
   dx = 1.0/(x(j)-x(j-1));
   dx1 = x(j)-x(j-1);
   dx2 = x(j)-x(j-2);
   u = V(:,j);
   % 1st order approximation
%   ux = dx*(V(:,jmax)-V(:,jmax-1));
   % 2nd order approximation
   ux = -(dx1*dx1*V(:,j-2)-dx2*dx2*V(:,j-1)-(dx1*dx1-dx2*dx2)*V(:,j)) ...
        /((dx1-dx2)*dx1*dx2);
   ut = Vt(:,j);
   Res(:,j) = feval(BC_R,t,xx,u,ux,ut);
         
   res = reshape(Res,nn,1);

% end of pdefdm_res

% --------------------------------------------------------------------------

function [SPDY,SPDYP] = jpattern(npde,jmax)
%
% this function defines the sparse pattern of the jacobian matrix.
%

   nn = npde*jmax;
   SPDY = sparse(nn,nn);
   SPDYP = sparse(nn,nn);
   
   for j=1:jmax
      j1 = j-1;
      j2 = j+1;
      if j==1, j1 = 1; j2 = 3; end
      if j==jmax, j1 = jmax-2; j2 = jmax; end
      SPDY((npde*(j-1)+1):(npde*j),(npde*(j1-1)+1):(npde*j2)) = ones(npde,(j2-j1+1)*npde);
   end
   for j=1:jmax
      SPDYP((npde*(j-1)+1):(npde*j),(npde*(j-1)+1):(npde*j)) = ones(npde,npde);
   end
   
% end of jpattern