% PDETWO - Solve two dimensional PDE's using ALGORITHM 565 PDETWO/PSETM/GEARB.
% 
% The types of PDE's that can be solved take the following form (in LaTeX notation)
%
% \begin{eqnarray*}
%  \frac{ \partial u_l }{ \partial t } &=& f_l( t,x,y,u_1,\ldots,u_{\mathrm{NPDE}},
%  \frac{ \partial u_1 }{ \partial x }, \ldots, \frac{ \partial u_{\mathrm{NPDE}} }{ \partial x },
%  \frac{ \partial u_1 }{ \partial y }, \ldots, \frac{ \partial u_{\mathrm{NPDE}} }{ \partial y }, \\
%  & & \frac{ \partial }{ \partial x } \left( \mathrm{DH}_{l,1} \frac{ \partial u_1 }{ \partial x } \right),\dots,\frac{ \partial }{ \partial x } \left( \mathrm{DH}_{l,\mathrm{NPDE}} \frac{ \partial u_{\mathrm{NPDE}} }{ \partial x } \right), \\ %\nonumber
%  & & \frac{ \partial }{ \partial y } \left( \mathrm{DV}_{l,1} \frac{ \partial u_1 }{ \partial y } \right),\dots,\frac{ \partial }{ \partial y } \left( \mathrm{DV}_{l,\mathrm{NPDE}} \frac{ \partial u_{\mathrm{NPDE}} }{ \partial y } \right) ) %\quad \mathrm{for} \quad l=1,2,\ldots,\mathrm{NPDE}
%  %\nonumber
% \end{eqnarray*}
% 
% With u a vector of NPDE unknown functions, u_x and u_y representing the 
% x and y spatical derivatives respectfully.  The "Y" boundary conditions 
% take the form of equations with the following form
% 
% \mathrm{AH}_l u_l + \mathrm{BH}_l \frac{ \partial u_l }{ \partial y } = \mathrm{CH}_l
% 
% while the "X" boundary conditions take the form of equations with the 
% following form
% 
% \mathrm{AV}_l u_l + \mathrm{BV}_l \frac{ \partial u_l }{ \partial x } = \mathrm{CV}_l
% 
% The boundary conditions need NOT be consistent with the initial conditions.  The 
% initial conditions are given by NPDE "functions" in the obvious form
% 
% u_l(t_0,x,y) = \phi_l(x,y) \quad \mathrm{for} \quad (x,y) \in R, \quad l=1,2,\ldots,\mathrm{NPDE}
% 
% The user of this routine must at specify the following octave functions:
% 
% A function to evaluate the right hand side f function (defined above)
% " """""""" "" """""""" """ horizontal diffusion coefficients DH(x,t,u)
% " """""""" "" """""""" """ vertical   diffusion coefficients DV(x,t,u)
% " """""""" "" """""""" """ horizontal boundary condition function AH(x,t,u), BH(x,t,u), and CH(x,t,u)
% " """""""" "" """""""" """ vertical   """""""" """"""""" """""""" AV(x,t,u), BV(x,t,u), and CV(x,t,u)
% A function to evaluate and return the initial conditions
% 
% The integration method is based on the method of lines (MOL).  The spatial 
% discretization is based on a finite difference discretization of the 
% two dimensional domain while the temporal integration is performed with the 
% GEARB ODE package.  The GEARB ODE package is nice in that it has methods for 
% integration both stiff and non-stiff problems.  Many more details can be
% found in the original article:
% 
% Algorithm 565: PDETWO/PSETM/GEARB: Solution of Systems of Two-Dimensional Nonlinear Partial Differential
% by David K. Melgaard, and Ricard F. Sincovec.
% ACM Transactions on Mathematical Software (TOMS), Volume 7 Issue 1 (March 1981)
% Pages: 126 - 135
% 
% Written by:
% -- 
% John L. Weatherwax                2006-03-05
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

close all; 

WL=-5.0; WR=+10.0; NX=2*30; 
H=10.0; NY=2*45; 

xgrid = linspace( WL, WR, NX ); ygrid = linspace( 0, H, NY ); 
tgrid = [ 0, 0.01 0.1:0.1:0.5 ];
U0 = UINIT( xgrid, ygrid );
[ U, opts ] = pdetwo("F","DIFF_H","DIFF_V","BNDRY_H","BNDRY_V", U0, ...
		     "xgrid", xgrid, "ygrid", ygrid, "tgrid", tgrid, "h", 1d-4, ...
		     "eps", 1d-4 );

tGrid = opts.tGrid; ntimes = length(tGrid); 
xbkpt = opts.xbkpt; 
ybkpt = opts.ybkpt;

[ X, Y ] = meshgrid(xbkpt,ybkpt);

% Only print the FINAL solution:
% 
for si=2:length(tgrid), % don't plot the initial condition ...
  u = squeeze(U(1,:,:,si)).';
  figure(si); mesh(X,Y,u); xlabel('X'); ylabel('Y'); zlabel('U'); title( 'approximate solution U' ); 
end

%save -mat eg1.mat tgrid X Y u;
save -mat eg1.mat tgrid xgrid ygrid X Y U;

return;