function acc05dim1fin
% One Dimensional Test for Distributed Parameter Numerics
% Some data from Wang 1999 JCR on BCNU study and Murray/Swanson99.
%%%%%
  clc   % clear variables;  caution: clears globals if after global command
%%%%%%
global RX RT RC
global alp12 kap13 alp21 kap23
global k1 k2 wy1 wy3 wu3
fprintf('\nacc05dim1fin.m:');
fprintf('\nOne Dimensional Test for Distributed Parameter Numerics');
np = 1; % Diagnostic print level: np=0 => none; np = 1 => some;
tolu = 0.5e-5; toly1 = 0.5e-7; toly2=0.5e-6;toly3=0.5e-3; % Shot Tolerances for U and Y;
tolcory = 0.5e-5; tolcorxi = 0.5e-5; % Correction Tolerances for Y and \xi;
Ncor = 5; % If big change, then use more corrections;
L = 10; % L double shots;
% (x,t) = (space,time); y(x,t) = (n1,n2,c) = (tumor,normal,c); xi(x,t) = co-state; 
DCT = 6.75e-5; DCN = 2.5e-6; % cm^2/sec; Wang99 BCNU Diffusion coefficients in Tumor/Normal
DTW = 4.2e-3; % cm^2/day; MurrayIIp567/Swanson99 Tumor diff. coef. in white matter; 
secsperday = 60^2*24; % Conversion seconds to days;
D = [DTW,1.0e-15,DCN*secsperday]; % Guess D(3) in Normal Tissue; 
% D(2)=nearzero due lack of data; Need D(1) for metastaces;
r = [0.1,0.0,0.0]; s = [0.0,0.0,0.2]; % Quadratic cost coefficients;
q = [0.1,0.0,0.1]; % FBH 9/08 Add quadratic final cost condition.
rho99 = 1.2e-2; % 1/day; Net growth rate by MurrayII/p567Swanson99, high grade glioma.
a3 = 1.31e-4; % (1/s): Effective absorption/reaction rate for BCNU, Wang99;
% Caution: Sign for a3 is in f3(y)=f(3,y);
a = [rho99,0.1e-10*secsperday,a3*secsperday]; % Guess normal growth rate a(2);
alp12 = 1.0e-4; kap13 = 1.0e-4; alp21 = 1.0e-4; kap23 = 1.0e-4; % All Guesses;
k1 = 1.0; k2 = 1.0; % FBH 9/08 Carrying Capacity Guesses as scale densities
wy1 = 1.e-3*k1; wy3 = 1.e-4; wu3 =1.e-5; %FBH 9/13; Gaussian Weights
T = 5; % Days
K = 1600; % 800; 400
dt = T/K; % Time step in days;
t = 0:dt:T;
RX = 1; % Domain Radius in cm;
RT = 2.e-2; % Guess resectioned tumor relative density or metastases extent;
RC = 1.e-2; % Guess catheter drug input;
M = [26]; % [14]; [010]; For explicit: keep K>10*max(D)*T*M^2/(2*RX)^2, 
%           so EffParMeshRatio=max(D)dt/dx^2<0.1, approximately;
dx = [2*RX/M(1)];
xt = -1:dx:1; x = xt';
u0 = u0f(x,t);
[mx,nx] = size(x); [mt,nt] = size(t); [lu0,mu0,nu0] = size(u0); 
ntm = round((nt+1)/2); ntq1 = round((ntm+1)/2); ntq3 = round((ntm+nt)/2);
mxm = round((mx+1)/2); mxq1 = round((mxm+1)/2); mxq3 = round((mxm+mx)/2);
fprintf('\ntolu=%6.1e; toly1=%6.1e; toly2=%6.1e; toly3=%6.1e;'...
    ,tolu,toly1,toly2,toly3);
fprintf('\nNote:  toly2 and toly3 not used.');
fprintf('\ntolcory=%6.1e; tolcorxi=%6.1e;',tolcory,tolcorxi);
fprintf('\nNcor=%i; L=%i',Ncor,L);
fprintf('\nM=%i; K=%i; mx=%i; nt=%i,mu0=%i; nu0=%i',M,K,mx,nt,mu0,nu0);
fprintf('\nT=%6.1e; RX(1)=%6.1e; RT=%6.1e; RC=%6.1e;',T,RX,RT,RC);
fprintf('\ndt=%6.1e; dx(1)=%6.1e; s(3)=%6.1e;',dt,dx(1),s(3));
fprintf('\nD=[%6.1e,%6.1e,%6.1e];',D);
fprintf('\na=[%6.1e,%6.1e,%6.1e];',a);
fprintf('\n[k1,k2]=[%6.1e,%6.1e];',k1,k2);
fprintf('\n[wy1,wy3,wu3]=[%6.1e,%6.1e,%6.1e];',wy1,wy3,wu3);
fprintf('\n[alp12,kap13,alp21,kap23]=[%6.1e,%6.1e,%6.1e,%6.1e];'...
    ,alp12,kap13,alp21,kap23);
fprintf('\nq=[%6.1e,%6.1e,%6.1e];',q);
fprintf('\nr=[%6.1e,%6.1e,%6.1e];',r);
fprintf('\ns=[%6.1e,%6.1e,%6.1e];',s);
fprintf('\nEPMR: max(D)*dt/dx(1)^2=%.1e;',max(D)*dt/dx(1)^2);
y0 = y0f(x); 
y(1:3,1:mx,1:nt) = zeros(3,mx,nt);
y(1:3,1:mx,1) = y0; % state IC
%%
if np >= 1 % Temp Print:
for i=1:3 
    for k = [1]
        fprintf('\nY(%i,x,%i)= ',i,k);
        for j = [1,2,3,4,mx]
            fprintf('%6.1e; ',y(i,j,k));
        end
    end
end, end
%%
for k = 1:nt
    xi(1:3,1:mx,k) = zeros(3,mx); % Co-State IC plus initial zeros for array;
    u(1:3,1:mx,k) = zeros(3,mx);
end % k
fprintf('\nmu3=[%i,%i,%i]; mxi3=[%i,%i,%i];\n'...
    ,size(u(3,1:mx,1:nt)),size(xi(3,1:mx,1:nt)));
ym(1:3,1:mx) = zeros(3,mx);
ddym(1:3,1:mx) = zeros(3,mx);
xim(1:3,1:mx) = zeros(3,mx);
ddxim(1:3,1:mx) = zeros(3,mx);
um(1:3,1:mx) = zeros(3,mx);
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Forward-Backward Shot Loop: Explicitized CNI first
Ncormaxy = 0; Ncormaxxi = 0; Lmax = 0;
U3old(1:mx,1:nt) = u0(3,1:mx,1:nt)+xi(3,1:mx,1:nt)/s(3);
Y1old(1:mx,1:nt) = y(1,1:mx,1:nt);
if np >= 1
    Y2old(1:mx,1:nt) = y(2,1:mx,1:nt);
    Y3old(1:mx,1:nt) = y(3,1:mx,1:nt);
end
for el = 1:L % Double Shot Loop
    % Backward Co-State part of Double Shot 
    % Regular Optimal Control from xi(3,*); %FH 9/10
    u(3,1:mx,1:nt) = U3old(1:mx,1:nt);
    y(1,1:mx,1:nt) = Y1old(1:mx,1:nt);
%%%%%%%%%%%%
if np >= 2 % Temp Print:    
for i=3 
    for k = [1,2,3,4,ntq1,ntm,ntq3,nt]
        fprintf('\nU(%i,x,%i)= ',i,k);
        for j = [1,2,3,4,mx]
            fprintf('%6.1e; ',u(i,j,k));
        end
    end
end, end
%%%%%%%%%%%%
    for k = 1:nt-1; 
        % Predictor (0th Corrector) Step:
        for j1 = 2:mx-1 % Get j=1&mx from BCs;
            ddy(1:3,j1) = (y(1:3,j1+1,k)-2*y(1:3,j1,k)+y(1:3,j1-1,k))/(dx(1))^2;
            Y = y(1:3,j1,k);
            for i = 1:3
                y(i,j1,k+1) = Y(i) + dt*(D(i)*ddy(i,j1)+a(i)*f(i,Y(i)) ...
                    + B(i,Y,t(k))*Y(i) + u(i,j1,k)); %FH 9/10
            end % i
        end % j1
        % No Flux BCs, not using D(i):
        y(1:3,mx,k+1) = (4*y(1:3,mx-1,k+1) - y(1:3,mx-2,k+1))/3;
        y(1:3,1,k+1) = (4*y(1:3,2,k+1) - y(1:3,3,k+1))/3;
        Y1kold(1:mx) = y(1,1:mx,k+1);
%%%%%%%%%%%%
if np >= 2 % Temp Print:
if k<11
    for i=1:3 
        for kp = [k+1]
            fprintf('\nY(%i,x,%i)= ',i,kp);
            for j = [1,2,3,4,mx]
                fprintf('%6.1e; ',y(i,j,kp));
            end
        end
    end
end, end
%%%%%%%%%%%%
        for ncor = 1:Ncor % Corrector Loop:
            % Predictor/Corrector Average Evaluation:
            ym(1:3,1:mx) = 0.5*(y(1:3,1:mx,k)+y(1:3,1:mx,k+1));
            tm = 0.5*(t(k)+t(k+1));
            for j1 = 2:mx-1
                ddym(1:3,j1) = (ym(1:3,j1+1)-2*ym(1:3,j1)+ym(1:3,j1-1))/(dx(1))^2;
            end % j1
            um(1:3,1:mx) = 0.5*(u(1:3,1:mx,k)+u(1:3,1:mx,k+1));
            % Corrector Step:
            for i = 1:3
                for j1 = 2:mx-1
                    Y = y(1:3,j1,k); Ym = ym(1:3,j1);
                    y(i,j1,k+1) = Y(i) + dt*(D(i)*ddym(i,j1)+a(i)*f(i,Ym(i)) ...
                        + B(i,Ym,tm)*Ym(i) + um(i,j1)); %FH 9/10
                end % j1
            end % i 
            % No Flux BCs, without D(i):
            y(1:3,mx,k+1) = (4*y(1:3,mx-1,k+1) - y(1:3,mx-2,k+1))/3;
            y(1:3,1,k+1) = (4*y(1:3,2,k+1) - y(1:3,3,k+1))/3;
%%%%%%%%%%%%
if np >= 2 % Temp Print:
    if k<11
        for i=1:3 
            for kp = [k+1]
                fprintf('\nY(%i,x,%i)= ',i,kp);
                for j = [1,2,3,4,mx]
                    fprintf('%6.1e; ',y(i,j,kp));
                end
            end
        end
    end 
end
%%%%%%%%%%%%
            Ncormaxy = max(Ncormaxy,ncor);
            Y1k(1:mx) = y(1,1:mx,k+1);
            if norm(Y1k(1:mx)-Y1kold(1:mx),inf) < tolcory*norm(Y1kold(1:mx),inf)
                break
            end
            Y1kold(1:mx) = y(1,1:mx,k+1);   
%%            
        end % ncor
    end % k and forward part of double shot
    % Backward Co-State part of Double Shot:
    % FBH 9/08
    xi(1,1:mx,nt) = -q(1)*y(1,1:mx,nt); % Co-State FC from terminal costs:
    xi(3,1:mx,nt) = -q(3)*y(3,1:mx,nt); % Co-State FC from terminal costs:
    %
    for k = nt:-1:2; 
        % Predictor (0th Corrector) Step:
        for j1 = 2:mx-1 % Get j=1&mx from BCs;
            ddxi(1:3,j1) = (xi(1:3,j1+1,k)-2*xi(1:3,j1,k)+xi(1:3,j1-1,k))/(dx(1))^2;
            XI = xi(1:3,j1,k); Y = y(1:3,j1,k); DB = dB(XI,Y,t(k));
            for i = 1:3
                xi(i,j1,k-1) = XI(i) + dt*(D(i)*ddxi(i,j1) ...
                    + df(i,Y(i))*a(i)*XI(i) - r(i)*Y(i) ...
                    + B(i,Y,t(k))*XI(i) + DB(i));
            end % i
        end % j1
        % No Flux BCs, without constant D(i):
        xi(1:3,mx,k-1) = (4*xi(1:3,mx-1,k-1) - xi(1:3,mx-2,k-1))/3;
        xi(1:3,1,k-1) = (4*xi(1:3,2,k-1) - xi(1:3,3,k-1))/3;
        XI1kold(1:mx) = xi(1,1:mx,k-1);
%%%%%%%%%%%%
if np >= 2 % Temp Print:
    if k>95
        for i=1:3 
            for kp = [k-1]
                fprintf('\nXI(%i,x,%i)= ',i,kp);
                for j = [1,2,3,4,mx]
                    fprintf('%6.1e; ',xi(i,j,kp));
                end
            end
        end
    end 
end
%%%%%%%%%%%%
        for ncor = 1:Ncor % Corrector Loop:
            % Predictor/Corrector Evaluation:
            xim(1:3,1:mx) = 0.5*(xi(1:3,1:mx,k)+xi(1:3,1:mx,k-1));
            ym(1:3,1:mx) = 0.5*(y(1:3,1:mx,k)+y(1:3,1:mx,k-1));
            for j1 = 2:mx-1
                ddxim(1:3,j1) = (xim(1:3,j1+1)-2*xim(1:3,j1)+xim(1:3,j1-1))/(dx(1))^2;
            end % j1
            tm = 0.5*(t(k)+t(k-1));
            % Corrector Step:
            for i = 1:3
                for j1 = 2:mx-1
                    XI = xi(1:3,j1,k); XIm = xim(1:3,j1); Ym = ym(1:3,j1);
                    DBm = dB(XIm,Ym,tm);
                    xi(i,j1,k-1) = XI(i) + dt*(D(i)*ddxim(i,j1) ...
                        + df(i,Ym(i))*a(i)*XIm(i) - r(i)*Ym(i) ...
                        + B(i,Ym,tm)*XIm(i) + DBm(i));
                end % j1
            end % i 
            % No Flux BCs, without D(i):
            xi(1:3,mx,k-1) = (4*xi(1:3,mx-1,k-1) - xi(1:3,mx-2,k-1))/3;
            xi(1:3,1,k-1) = (4*xi(1:3,2,k-1) - xi(1:3,3,k-1))/3;
%%%%%%%%%%%%
if np >= 2 % Temp Print:
    if k>95
        for i=1:3 
            for kp = [k-1]
                fprintf('\nXI(%i,x,%i)= ',i,kp);
                for j = [1,2,3,4,mx]
                    fprintf('%6.1e; ',xi(i,j,kp));
                end
            end
        end
    end 
end
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
            Ncormaxxi = max(Ncormaxxi,ncor);
            XI1k(1:mx) = xi(1,1:mx,k-1);
            if norm(XI1k(1:mx)-XI1kold(1:mx),inf) < tolcorxi*norm(XI1kold(1:mx),inf)
                break
            end
            XI1kold(1:mx) = xi(1,1:mx,k-1);
        end % ncor
    end % k and backward part of double shot 
    Lmax = max(Lmax,el);
    % Final Regular Optimal Control from xi(3,*);
    u(3,1:mx,1:nt) = u0(3,1:mx,1:nt)+xi(3,1:mx,1:nt)/s(3);
    U3(1:mx,1:nt) = u(3,1:mx,1:nt); 
    Y1(1:mx,1:nt) = y(1,1:mx,1:nt);
    U3oldnorm = norm(U3old(1:mx,1:nt),inf);
    DU3norm = norm(U3(1:mx,1:nt)-U3old(1:mx,1:nt),inf);
    Y1oldnorm = norm(Y1old(1:mx,1:nt),inf);
    DY1norm = norm(Y1(1:mx,1:nt)-Y1old(1:mx,1:nt),inf);
    U3old(1:mx,1:nt) = U3(1:mx,1:nt);
    Y1old(1:mx,1:nt) = Y1(1:mx,1:nt);
    if np >= 1
        Y2(1:mx,1:nt) = y(2,1:mx,1:nt);
        Y3(1:mx,1:nt) = y(3,1:mx,1:nt);
        Y2oldnorm = norm(Y2old(1:mx,1:nt),inf);
        DY2norm = norm(Y2(1:mx,1:nt)-Y2old(1:mx,1:nt),inf);
        Y3oldnorm = norm(Y3old(1:mx,1:nt),inf);
        DY3norm = norm(Y3(1:mx,1:nt)-Y3old(1:mx,1:nt),inf);
        Y2old(1:mx,1:nt) = Y2(1:mx,1:nt);
        Y3old(1:mx,1:nt) = Y3(1:mx,1:nt);
    end
    if (DU3norm < tolu*U3oldnorm) & (DY1norm < toly1*Y1oldnorm)
        break
    end % if
end % el shot loop
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%
fprintf('\nNcormaxy=%i; Ncormaxxi=%i; Lmax=%i;',Ncormaxy,Ncormaxxi,Lmax);
fprintf('\nDU3norm=%6.1e; tolu*U3oldnorm=%6.1e; U3oldnorm=%6.1e;'...
    ,DU3norm,tolu*U3oldnorm,U3oldnorm);
fprintf('\nDY1norm=%6.1e; toly1*Y1oldnorm=%6.1e; Y1oldnorm=%6.1e;\n'...
        ,DY1norm,toly1*Y1oldnorm,Y1oldnorm);
%%
if np >= 1
    fprintf('\nDY2norm=%6.1e; toly2*Y2oldnorm=%6.1e; Y2oldnorm=%6.1e;'...
        ,DY2norm,toly2*Y2oldnorm,Y2oldnorm);
    fprintf('\nDY3norm=%6.1e; toly3*Y3oldnorm=%6.1e; Y3oldnorm=%6.1e;\n'...
        ,DY3norm,toly3*Y3oldnorm,Y3oldnorm);
end
%%
if np >= 2 % Temp Print:
for i=3 
    for k = [1,2,3,4,ntq1,ntm,ntq3,nt]
        fprintf('\nU(%i,x,%i)= ',i,k);
        for j = [1,2,3,mxq1,mxm,mxq3,mx]
            fprintf('%6.1e; ',u(i,j,k));
        end
    end
end, end
% Temp Print:
if np >= 1 % Temp Print:
for i=1:1 %for i=1:3
    for k = [1,2,3,4,ntq1,ntm,ntq3,nt]
        fprintf('\nY(%i,x,%i)= ',i,k);
        for j = [1,2,3,mxq1,mxm,mxq3,mx]
            fprintf('%6.1e; ',y(i,j,k));
        end
    end
    fprintf('\n100*y(1,mxm,nt)/y(1,mxm,1) = %6.2f%%;\n',100*y(1,mxm,nt)/y(1,mxm,1));
end,end
% Temp Print:
if np >= 2 % Temp Print:
for i=1:3 
    for k = [nt,nt-1,nt-2,ntq3,ntm,ntq1,1]
        fprintf('\nxi(%i,x,%i)= ',i,k);
        for j = [1,2,3,mxq1,mxm,mxq3,mx]
            fprintf('%6.1e; ',xi(i,j,k));
        end
    end
end, end
fprintf('\n');
%%
nfig=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
nfig = nfig +1;
scrsize = get(0,'ScreenSize');
ss = [4.0,3.8,3.6,3.4,3.2];
figure(nfig) % 
plot(x,u(3,1:mx,1),'k:o',x,u(3,1:mx,ntq1),'k:',x,u(3,1:mx,ntm),'k-.' ...
    ,x,u(3,1:mx,ntq3),'k--',x,u(3,1:mx,nt),'k-','linewidth',4);
axis tight
%grid on
hlegend = legend('U_3^*(x,0)','U_3^*(x,t_{q1})','U_3^*(x,t_{mid})'...
    ,'U_3^*(x,t_{q3})','U_3^*(x,t_f)',0);
htitle = title('Optimal Control U_3^*(x,t)');
hxlabel = xlabel('x, Space in cm.');
hylabel = ylabel('U_3^*(x,t), Regular Control');
set(hlegend,'Fontsize',32,'FontName','Helvetica','FontWeight','Bold');
set(gca,'Fontsize',28,'FontName','Helvetica','FontWeight','Bold','linewidth',3);
set(htitle,'Fontsize',36,'FontName','Helvetica','FontWeight','Bold');
set(hxlabel,'Fontsize',32,'FontName','Helvetica','FontWeight','Bold');
set(hylabel,'Fontsize',32,'FontName','Helvetica','FontWeight','Bold');
set(gcf,'Color','White','Position' ...
    ,[scrsize(3)/ss(nfig) 35 scrsize(3)*0.70 scrsize(4)*0.86]);
%%%%%%%%%%%%%%%%%%%%
nfig = nfig +1;
figure(nfig) % 
% plot(x,y(3,1:mx,1),'k:o',...
plot(...
    x,y(3,1:mx,ntq1),'k:',x,y(3,1:mx,ntm),'k-.' ...
    ,x,y(3,1:mx,ntq3),'k--',x,y(3,1:mx,nt),'k-','linewidth',5);
axis tight
%grid on
% hlegend = legend('C^*(x,0)',...
hlegend = legend(...
    'C^*(x,t_{q1})','C^*(x,t_{mid})','C^*(x,t_{q3})','C^*(x,t_f)',0);
htitle = title('Optimal Drug Concentration, C^*(x,t)');
hxlabel = xlabel('x, Space in cm.');
hylabel = ylabel('C^*(x,t), Drug Concentration');
set(hlegend,'Fontsize',32,'FontName','Helvetica','FontWeight','Bold');
set(gca,'Fontsize',28,'FontName','Helvetica','FontWeight','Bold','linewidth',3);
set(htitle,'Fontsize',36,'FontName','Helvetica','FontWeight','Bold');
set(hxlabel,'Fontsize',32,'FontName','Helvetica','FontWeight','Bold');
set(hylabel,'Fontsize',32,'FontName','Helvetica','FontWeight','Bold');
set(gcf,'Color','White','Position' ...
    ,[scrsize(3)/ss(nfig) 35 scrsize(3)*0.70 scrsize(4)*0.86]);
%%%
%%%%%%%%%%%%%%%%%%%%
nfig = nfig +1;
figure(nfig) % 
plot(x,y(1,1:mx,1),'k:o',x,y(1,1:mx,ntq1),'k:',x,y(1,1:mx,ntm),'k-.' ...
    ,x,y(1,1:mx,ntq3),'k--',x,y(1,1:mx,nt),'k-','linewidth',5);
axis tight
%grid on
hlegend = legend('N_1^*(x,0)','N_1^*(x,t_{q1})'...
    ,'N_1^*(x,t_{mid})','N_1^*(x,t_{q3})' ...
    ,'N_1^*(x,t_f)',0);
htitle = title('Optimal Relative Tumor Density N_1^*(x,t)');
hxlabel = xlabel('x, Space in cm.');
hylabel = ylabel('N_1^*(x,t)');
set(hlegend,'Fontsize',32,'FontName','Helvetica','FontWeight','Bold');
set(gca,'Fontsize',28,'FontName','Helvetica','FontWeight','Bold','linewidth',3);
set(htitle,'Fontsize',36,'FontName','Helvetica','FontWeight','Bold');
set(hxlabel,'Fontsize',32,'FontName','Helvetica','FontWeight','Bold');
set(hylabel,'Fontsize',32,'FontName','Helvetica','FontWeight','Bold');
set(gcf,'Color','White','Position' ...
    ,[scrsize(3)/ss(nfig) 35 scrsize(3)*0.70 scrsize(4)*0.86]);
%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Subfunctions:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y0v = y0f(x)
% initial plus boundary states:
global RX RT RC
global k1 k2 wy1 wy3 wu3
m01 = 0.0; s01 = RT;
m03 = 0.0; s03 = RC;
[mx,nx] = size(x); xones = ones(size(x))';% FBH 9/08
y0v(1,1:mx) = 0.0*xones; % Pre-Initialization for Tumor size;
y0v(2,1:mx) = k1*xones; % Scale population Density to Normal Density;
y0v(3,1:mx) = 0.0*xones; % Guess zero drug start;
for j = 2:mx-1 % FBH 9/09 Gaussian IC after MurrayII; Detectable 4e+4 cells/cm^2?
    y0v(1,j) = wy1*exp(-0.5*((x(j,1)-m01)/s01)^2)/sqrt(2*pi*s01^2);
    y0v(3,j) = wy3*exp(-0.5*((x(j,1)-m03)/s03)^2)/sqrt(2*pi*s03^2);
%    if abs(x(i,1)) < RT
%        y0v(1,i) = 1.e-8*k1; % FBH 9/08 Guess relative initial Tumor cell Density;
%    end
end
% No Flux Boundary correction compatibility:
y0v(1:3,mx) = (4*y0v(1:3,mx-1)-y0v(1:3,mx-2))/3;
y0v(1:3,1) = (4*y0v(1:3,2)-y0v(1:3,3))/3; 
%% End y0f
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function u0v = u0f(x,t) % FBH 9/10 new:
% Threshold drug control:
global RX RT RC
global k1 k2 wy1 wy3 wu3
m03 = 0.0; s03 = RC;
[mx,nx] = size(x);% FBH 9/08
[mt,nt] = size(t);
for k = 1:nt
    for j = 2:mx-1 % FBH 9/09 Gaussian IC after MurrayII; Detectable 4e+4 cells/cm^2?
        u0v(3,j,k) = wu3*exp(-0.5*((x(j,1)-m03)/s03)^2)/sqrt(2*pi*s03^2);
    end
% No Flux Boundary correction compatibility:
    u0v(3,mx,k) = (4*u0v(3,mx-1,k)-u0v(3,mx-2,k))/3;
    u0v(3,1,k) = (4*u0v(3,2,k)-u0v(3,3,k))/3; 
end
u0v(1:2,1:mx,1:nt) = 0.0;
%% End u0f
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function fv = f(i,yi)
% Growth functions, less net growth rates (assume diagonal models):
global k1 k2 wy1 wy3 wu3
if i==1
% Tumor growth: assume exponential -> logistic:
    fv = + yi*(1-yi/k1);
elseif i==2
% Normal Growth: Assume logistic rate, since no carrying capacity K2 available;
    fv = + yi*(1-yi/k2);
elseif i==3
% Drug absorption/reaction: assume negative exponential:
    fv = - yi; % Note minus sign here, and not on a3;
end
%% End f
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dfv = df(i,yi)
global k1 k2 wy1 wy3 wu3
% Growth function gradient, less net growth rates (assume diagonal models: & j==i):
if i==1 
% Tumor growth: assume exponential
    dfv = + 1-2*yi/k1;
elseif i==2
% Normal Growth: Assume exponential rate, since no carrying capacity K2 available;
    dfv = + 1-2*yi/k2;
elseif i==3
% Drug absorption/reaction: assume negative exponential:
    dfv = - 1;
end
%% End df
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function Bv = B(i,Y,t)
% Nonlinear Coefficient functions (assume diagonal models):
global alp12 kap13 alp21 kap23
if i==1
% Tumor growth: assume exponential:
    Bv = -(alp12*Y(2)+kap13*Y(3));
elseif i==2
% Normal Growth: Assume exponential rate, since no carrying capacity K2 available;
    Bv = -(alp21*Y(1)+kap23*Y(3));
elseif i==3
% Drug absorption/reaction: assume negative exponential:
    Bv = 0;
end
%% End B
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dBv = dB(XI,Y,t)
% Nonlinear Coefficient Gradient grad_{Y_i)[B_j](Y,t)*XI_j*Y_j (assume diagonal models):
global alp12 kap13 alp21 kap23
%
dBv(1) = -alp21*XI(2)*Y(2);
%
dBv(2) = -alp12*XI(1)*Y(1);
%
dBv(3) = -(kap13*XI(1)*Y(1)+kap23*XI(2)*Y(2));
% 
%% End dB
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% End acc05dim1sf.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%