function ppde2dbc
% ppde-dbc.m MCS571 Test  for u_t= A*u_{xx}, u_x(a,t)=0; 
%        u_x(b,t)=-B*(u-TA)/K; u(x,0)=0
% 
clc % clear variables, but must come before globals, 
clf % clear figures
%
fprintf('\nDerivative BCs PPDE by 2nd Order FFD/BFD Treatment (MATLAB6.5):')
scrsize = get(0,'ScreenSize');
ss = [3.3,3.1,2.9];
f=0;
m = 7; n=15; a=0; b=2; h = (b-a)/(m+1); A=1; B=1; K=1; BDK=B/K; TA=100; 
Den = 3+2*h*BDK; Num=(1+2*h*BDK); 
X=h*(0:m+1)';
U=zeros(m+2,n+1); % Global U with (1,1) subscript base in stead of (0,0)
U(:,1) = u0(X); j=0;
U(1,1) = 0.5*((4*U(2,1)-U(3,1))/3+U(1,1)); %Average Condition at corner
U(m+2,1) = 0.5*((4*U(m+1,1)-U(m,1)+2*h*BDK*TA)/Den+U(m+2,1)); 
% Average Condition at corner
fprintf('\nj=%2i: [%f %f %f %f %f %f %f %f %f]',j,U(:,j+1))
rv=[0.25,0.50,0.65]; % Marginally Stable, Stable,  Unstable Cases;
stitle = {'PPDE EEM: u_t=A*u_{xx}, r = 0.25, Stable' ...
        ,'PPDE EEM: u_t=A*u_{xx}, r = 0.5, Still Stable' ...
        ,'PPDE EEM: u_t=A*u_{xx}, r = 0.65, Unstable'};
for ir = 1:3
    r=rv(ir); 
    k=r*h^2; ar = A*r; ar2=ar*2; ar2m = ar2*Num/Den;
fprintf('\n  m=%i; n=%i; h=%f; r=%f; k=%f; 2*A*r=%f; 2*A*r*Num/Den=%f;' ...
    ,m,n,h,r,k,ar2,ar2m)
    for j=1:n
        U(2:m+1,j+1) = (1-2*ar)*U(2:m+1,j)+ar*(U(3:m+2,j)+U(1:m,j)); 
%       Usual Euler FFD(t)=CFD(x)
        U(1,j+1) = (4*U(2,j+1)-U(3,j+1))/3; % FFD^(2) LBC
        U(m+2,j+1) = (4*U(m+1,j+1)-U(m,j+1)+2*h*BDK*TA)/Den;  %BFD^(2) RBC
        fprintf('\nj=%2i: [%f %f %f %f %f %f %f %f %f]',j,U(:,j+1))
    end
    f=f+1;
figure(f)
plot(X,U(:,1),'b-',X,U(:,2),'b-',X,U(:,3),'b-',X,U(:,4),'b-',X,U(:,5),'b-' ...
    ,X,U(:,6),'b-',X,U(:,7),'b-',X,U(:,8),'b-',X,U(:,9),'b-',X,U(:,10),'b-' ...
    ,X,U(:,11),'b-',X,U(:,12),'b-',X,U(:,13),'b-',X,U(:,14),'b-',X,U(:,15),'b-' ...
    ,X,U(:,16),'b-','LineWidth',2) % max input arguments!
    axis tight
    htitle=title(stitle(ir));
    hxlabel = xlabel('x_i=a+i*h,  i=0:m+1');
    hylabel = ylabel('U(x_i,t_j),  j=0:n');
    htext = text(0.25,37,'Derivative BC by 2nd Order FFD/BFD');
    set(gca,'Fontsize',20,'FontName','Helvetica','FontWeight','Bold' ...
        ,'linewidth',2);
set(htitle,'Fontsize',20,'FontName','Times New Roman','FontWeight','Bold')
set(hxlabel,'Fontsize',20,'FontName','Times New Roman','FontWeight','Bold')
set(hylabel,'Fontsize',20,'FontName','Times New Roman','FontWeight','Bold')
set(htext,'Fontsize',16,'FontName','Times New Roman','FontWeight','Bold')
    set(gcf,'Color','White','Position' ...
        ,[scrsize(3)/ss(ir) 35 scrsize(3)*0.60 scrsize(4)*0.86]);
%    MATLAB defines the figure Position property as a vector. 
%            [left bottom width height]
end % r for loop
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function u0v = u0(x)
u0v = 0;
% END u0