% PDECOL - Solve one dimensional PDE's using ALGORITHM 540 PDECOL.
% 
% The types of PDE's that can be solved take the following form
% 
% u_t = f(t,x,u,u_x,u_xx)
% 
% With u a vector of NPDE unknown functions, u_t and u_x representing time
% and spacial derivatives respectfully.  Boundary conditions take the form
% of equations of the following form
% 
% b(u,u_x) = z(t)
% 
% and must be consistent with the initial conditions.  The user of the 
% routine must at minimum specify the following octave functions:
% 
% A function to evaluate f (defined above)
% A function to evaluate the derivative of b with respect to u
% " """""""" "" """""""" """ """""""""" "" " """" """"""" "" u_x
% " """""""" "" """""""" """ """""""""" "" z """" """"""" "" t
% A function to evaluate and return the initial conditions
% 
% This method in space is a finite element collocation method with piecewise 
% polynomials (B-Splines).  The time integration is performed using the method
% of lines with choices for several standard ordinary differential equations 
% solvers.  Many more details can be found in the original article:
% 
% Algorithm 540: PDECOL, General Collocation Software for Partial Differential Equations
% by N. K. Madsen, R. F. Sincovec
% ACM Transactions on Mathematical Software (TOMS), Volume 5 Issue 3 (September 1979)
% 
% Written by:
% -- 
% John L. Weatherwax                2005-04-27
% 
% email: wax@alum.mit.edu
% 
% Please send comments and especially bug reports to the
% above email address.
% 
%-----

if( 0 ) 
  [ U, opts ] = pdecol("F","DBDU_BNDRY","DBDUX_BNDRY","DZDT_BNDRY","UINIT" );
else
  [ U, opts ] = pdecol( "F","DBDU_BNDRY","DBDUX_BNDRY","DZDT_BNDRY","UINIT",...
			"dfdu", "DFDU_DERIVF","dfdux", "DFDUX_DERIVF", "dfduxx", "DFDUXX_DERIVF", ...
			"MF", 21, "tgrid", linspace(0,0.1,50), 'xgrid', linspace(0,1,100) );
end

tGrid = opts.tGrid; 
xbkpt = opts.xbkpt; 

%tGrid(1) = []; 
[ X, T ] = meshgrid(xbkpt,tGrid);
u = squeeze(U(1,:,:)).';
v = squeeze(U(2,:,:)).';
figure; mesh(X,T,u); xlabel('X'); ylabel('T'); zlabel('U'); 
figure; mesh(X,T,v); xlabel('X'); ylabel('T'); zlabel('V'); 

return; 

close all; 
ntouts = length(tGrid)-1; 
h=figure; 
for nti=1:ntouts,
  figure(h); 
  subplot( 1,2,1 ); %axis( [ 0 1 0.4 1.8 ] ); 
  plot( xbkpt, U(1,:,nti) ); grid; %title( [ 'U(1) at t(',num2str(nti),')' ] ); 
  %pause; 
  subplot( 1,2,2 ); %axis( [ 0 1 2.7 3.7 ] ); 
  plot( xbkpt, U(2,:,nti) ); grid; %title( [ 'U(2) at t(',num2str(nti),')' ] ); 
  pause;
end