%
% Dynamic analysis of drug delivery
%
% The following equations model a polymer matrix in cylindrical
% coordinates 
%
%   u1 material balance
%
%   u1t = Du1*(u1rr + (1/r)u1r)                                 (1)
%
%   u2 material balance
%
%   u2t = Du2*(u2rr + (1/r)*u2r)                                (2)
%
% The variables and parameters for this model are (in cgs units)
%
%   Water concentration            u1                     (eq. (1))
%
%   Drug concentration             u2                     (eq. (2))
%
%   Time                            t
%
%   Radial position                 r
%
%   Axial position                  z
%
%   Polymer matrix radius          r0              1
%
%   Polymer matrix length          zL              2
%
%   Initial u1                    u10            0.5
%
%   Initial u2                    u20              1
%
%   External u1                   u1e              1
%
%   External u2                   u2e              0
%
%   u1 diffusivity                Du1        1.0e-06
%
%   u2 diffusivity                Du2        1.0e-06
%
%   u1 mass transfer coefficient  ku1            0.1
%
%   u2 mass transfer coefficient  ku2            0.1
%
% The solution to this system, u1(r,z,t) (from eq. (1)) 
% and u2(r,z,t) (from eq. (2)) is computed by ode15s for 
% the integration in t. 
%
% Clear previous files
  clear all
  clc
%
% Global area
%
% Global area
  global     nr     nz     dr     dz    drs    dzs...
              r     r0      z     zL    Du1    Du2...
             u1     u2    u1e    u2e    ku1    ku2...
          ncall
%
% Model parameters
  u10=0.5; u20=1;
  u1e=1; u2e=0;
  r0=1; zL=2;
  Du1=1.0e-06; Du2=1.0e-06;
  ku1=1.0e-01; ku2=1.0e-01;
%
% Grid in axial direction
  nz=11;
  dz=zL/(2*(nz-1));
  for i=1:nz
    z(i)=zL/2+(i-1)*dz;
  end
  dzs=dz^2;
%
% Grid in radial direction
  nr=11;
  dr=r0/(nr-1);
  for j=1:nr
    r(j)=(j-1)*dr;
  end
  drs=dr^2;  
%
% Independent variable for ODE integration
  tf=2*3600*24;
  tout=[0.0:2*900*24:tf]'; 
  nout=5;
  ncall=0;
%
% Initial condition
  for i=1:nz
  for j=1:nr
    u1(i,j)=u10;
    u2(i,j)=u20;
    u0((i-1)*nr+j)=u1(i,j);
    u0((i-1)*nr+j+nz*nr)=u2(i,j);
  end
  end
%
% ODE integration 
  reltol=1.0e-04; abstol=1.0e-04;
  options=odeset('RelTol',reltol,'AbsTol',abstol);
  mf=1;
  if(mf==1)[t,u]=ode15s(@pde_1,tout,u0,options);end 
  if(mf==2)[t,u]=ode15s(@pde_2,tout,u0,options);end 
  if(mf==3)[t,u]=ode15s(@pde_3,tout,u0,options);end  
  if(mf==4)[t,u]=ode15s(@pde_4,tout,u0,options);end   
%
% 1D to 2D matrices
  for it=1:nout
  for i=1:nz
  for j=1:nr  
    u1(it,i,j)=u(it,(i-1)*nr+j);
    u2(it,i,j)=u(it,(i-1)*nr+j+nz*nr);
  end
  end
  end
%
% Display a heading and radial profiles
%
% z=zL/2
  fprintf('\n  nr = %2d   nz = %2d\n',nr,nz);
  for it=1:nout
    fprintf('\n t = %4.1f   z = %3.1f\n',t(it)/3600,z(1));  
  for j=1:nr
    fprintf(...
        ' r = %4.2f   u1(r,z,t) = %6.3f   u2(r,z,t) = %6.3f\n',...
        r(j),u1(it,1,j),u2(it,1,j));
  end
    fprintf('\n ku2*(u2(r=r0,z=zL/2,t)-u2e) = %8.3e\n\n',...
            3600*ku2*(u2(it,1,nr)-u2e));
  end
  fprintf('\n ncall = %5d\n',ncall);
%
% Parametric plots, radial profiles, z=zL/2
%
% 2D plots
  figure(1);
  subplot(2,2,1)
  for j=1:nr u1plot(j)=u1(2,1,j); u2plot(j)=u2(2,1,j); end
  plot(r,u1plot,'-o',r,u2plot,'-x'); axis([0 r0 0 1]);
  title('u1,u2(r,z=zL/2,t=12hr),o-u1,x-u2');
  xlabel('r'); ylabel('u1,u2(r,z=zL/2,t=12hr)')
  subplot(2,2,2)
  for j=1:nr u1plot(j)=u1(3,1,j); u2plot(j)=u2(3,1,j); end
  plot(r,u1plot,'-o',r,u2plot,'-x'); axis([0 r0 0 1]);
  title('u1,u2(r,z=zL/2,t=24hr),o-u1,x-u2');  
  xlabel('r'); ylabel('u1,u2(r,z=zL/2,t=24hr)')
  subplot(2,2,3)
  for j=1:nr u1plot(j)=u1(4,1,j); u2plot(j)=u2(4,1,j); end
  plot(r,u1plot,'-o',r,u2plot,'-x'); axis([0 r0 0 1]);
  title('u1,u2(r,z=zL/2,t=36hr),o-u1,x-u2');   
  xlabel('r'); ylabel('u1,u2(r,z=zL/2,t=36hr)')
  subplot(2,2,4)
  for j=1:nr u1plot(j)=u1(5,1,j); u2plot(j)=u2(5,1,j); end
  plot(r,u1plot,'-o',r,u2plot,'-x'); axis([0 r0 0 1]);
  title('u1,u2(r,z=zL/2,t=48hr),o-u1,x-u2');   
  xlabel('r'); ylabel('u1,u2(r,z=zL/2,t=48hr)')
%
% 3D plots
  for it=1:nout
    for j=1:nr u1surf(it,j)=u1(it,1,j); end
    for j=1:nr u2surf(it,j)=u2(it,1,j); end    
  end
  figure(2)
  surf(r,t/(3600*24),u1surf);
  xlabel('r (cm)'); ylabel('t (day)'); zlabel('u1(r,z=zL/2,t)');
  figure(3)
  surf(r,t/(3600*24),u2surf);
  view(-20,60)
  xlabel('r (cm)'); ylabel('t (day)'); zlabel('u2(r,z=zL/2,t)');  
