# Supplementary Material (ESI) for Physical Chemistry Chemical Physics
# This journal is (c) The Owner Societies 2011

%fullmodel081110_V2.m

%--------------------------------------------------------------------------
% Finite difference solver for diffusion-migration with a conducting polymer
% electrode.
% Version 2.0
% Karthik Kannappan           Gib Bogle
% k.kannappan@auckland.ac.nz  g.bogle@auckland.ac.nz
%--------------------------------------------------------------------------
% The smoothing function, fastsmooth.m, can be downloaded from:
% http://www.mathworks.com/matlabcentral/fileexchange/19998-fast-smoothing-function
%--------------------------------------------------------------------------

function ep1
clear all
% format long e

% Constants
RC=8.314472;                    % Gas constant ( J K-1 mol-1 )
TC=300;                         % Room Temperature ( K )
FC=96485.3383;                  % Faraday's constant  ( C mol-1 )
charge= 1;                      % Charge of the species
DA= 2.0e-9;                     % Diffusion constant A ( m2 s-1 ) (2.0e-9)
DB= DA;                         % Diffusion constant B ( m2 s-1 )
UA= FC*DA/(RC*TC);              % Mobility of Ion A ( m2 s-1 v-1 )
UB= FC*DB/(RC*TC);              % Mobility of Ion B  ( m2 s-1 v-1 )

C_F_L = 0.005;                  % Courant-Friedrichs-Levy constant (actually C/u)

% Non-dimensional constants
C0 = 10;                        % Bulk Concentration (mol / m3)= 0.01 mol/liter
V0= RC*TC/FC;                   % 0.025852 Volts 
X0= 0.001;                      % Distance between electrodes (m) (0.001)
T0= (X0*X0)/DA;                 % Time scale of diffusion
R0 = 1/((UA+UB)*C0*FC*X0)       % reference resistance (ohm)
I0 = (UA+UB)*C0*RC*TC/X0;       % reference value for current density (IS)
i0 = I0*X0*X0;                  % the reference value to non-dimensionalise IP is I0*X0*X0

% Geometry of the cell
XDistance= 10;                  % Non-dimensional
YDistance= 1;                   % Non-dimensional

%Grid
nx=100;
ny=30; 
N=nx*ny;
dx=XDistance/(nx-1);
dy=YDistance/(ny-1);

%Length of the counter-electrode 
Length=1;
Variable=nx-(fix(nx*Length));
OXinitial = 0.01;

% Parameter values (actual units)
Mion_actual = 1.0;                                    % Ion capacity of the film % (moles / m3)
STDP_actual = V0*log((1-OXinitial)/OXinitial);        % Standard potential (V)
Width_actual = 0.001;                                 % Width of the polymer (m)
Thickness_actual = 0.00001;                           % Thickness of the polymer (m)

MT = Mion_actual*Thickness_actual;

emin = OXinitial;
emax = 0.995;

% Change to initial applied potential set to generate zero current
Vinitial = STDP_actual - V0*log((1-emin)/emin);
Vfinal = STDP_actual - V0*log((1-emax)/emax);

Nramp = 20;
dVdt = (Vfinal - Vinitial)/Nramp;

dt = C_F_L/((1/dx) + (1/dy));

simtime = input('Enter simulation time in seconds ');
nt = simtime/(dt*T0) + 1;

nplot=fix(nt/10);
 
% Applied parameters ( non dimensionalised )
STDP = STDP_actual/V0;              % Standard potential                                         
Width = Width_actual/X0;            % Width of the polymer 
Thickness = Thickness_actual/X0;    % Thickness of the polymer 
AreaP = Width*Thickness;            % Area of the polymer cross section 
AreaS = Width*dx;                   % Area of polymer-solution interface 

for ix=1:nx    
    Bpottop(ix) =0.0/V0;     % Counter electrode Potential 
    F(ix)= 0;                % Flux at the boundaries
end

% Parameters defining electrolyte properties

%%%%%%%%%%%%%%%%%%%
%Initial conditions
%%%%%%%%%%%%%%%%%%%

OX(1:nx) = OXinitial;

%Initial concentration of the solution
for ix = 1:nx
    for iy = 1:ny
        x(ix)= (ix-1)*dx;
        y(iy)= (iy-1)*dy;
        C(ix,iy) = 1;       % Non-dimensional   
    end
end

ESigma=(UA+UB)*FC*C0   %( ohm-1 m-1)
Sigmaratio= input('Enter the electrolyte to polymer conductivity ratio');
% Parameters defining polymer properties
Sigmamin=ESigma/Sigmaratio    % Minimum conductivity of the polymer  (S/m)
Sigmamax=Sigmamin*10000       % Maximum conductivity of the polymer  (S/m)
OXmax=1;                      % Maximum oxidation state
OXT=0.5*OXmax;                % Oxidation state at which transition occurs
w=0.03*OXmax;                 % Width of transition of polymer from insulating to conducting state

Sigmaratio=ESigma/Sigmamin

% initial condition display
figure(1)
 mesh(C)
 xlabel('y')
 ylabel('x')
 title('Initial condition for C')
 pause;
  
% Starting time iteration

for t = 1:nt

% ramp applied voltage
    if (dVdt > 0)
        PPB_actual = min(Vinitial + dVdt*t,Vfinal);
    else
        PPB_actual = max(Vinitial + dVdt*t,Vfinal);
    end
    PPB = PPB_actual/V0;                                                    % applied potential (non-dimensional)
       
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Part 1 : Finding polymer resistivity and delta potential
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    for ix=1:nx
        Sigma(ix)=(Sigmamax-Sigmamin)/(1+exp(-(OX(ix)-OXT)/w))+Sigmamin;   % Actual conductivity at the mesh points ( ohm-1 m-1)
    end
    for ix=2:nx
        RPOLYMER(ix)=((1/Sigma(ix) + 1/Sigma(ix-1))/2)*((dx*X0)/(Thickness_actual*Width_actual)); % Actual average resistivity b/w two mesh points ( ohms)
    end
    RPOLYMER(1) = 0;
    
    RP = RPOLYMER/R0;                                                      % Non-dimensional average resistivity

    for ix=1:nx
         Pdel(ix)=STDP - log(C(ix,1)*(1-OX(ix))/OX(ix));                   % Delta potential (non-dimensional )
    end
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % End of Part 1
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Part 2 : Finding potential gradient in the solution
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%Concentration derivatives%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % derivatives of the concentration at each points
    for ix = 2:nx-1
        for iy = 2:ny-1
            cfx(ix,iy)=(C(ix+1,iy)-C(ix-1,iy))/(2*dx);
            cfy(ix,iy)=(C(ix,iy+1)-C(ix,iy-1))/(2*dy);
            csxy(ix,iy)=(C(ix+1,iy)+C(ix-1,iy)-2*C(ix,iy))/(dx*dx)+ (C(ix,iy+1)+C(ix,iy-1)-2*C(ix,iy))/(dy*dy);
        end
    end
    %left and right boundary
    for iy = 2:ny-1
        %first derivative of ix
        cfx(1,iy)=(C(2,iy)-C(1,iy))/(dx); 
        cfx(nx,iy)=(C(nx,iy)-C(nx-1,iy))/(dx);
        %first derivative of iy
        cfy(1,iy)=(C(1,iy+1)-C(1,iy-1))/(2*dy);
        cfy(nx,iy)=(C(nx,iy+1)-C(nx,iy-1))/(2*dy);     
        %second derivative x,y
        csxy(1,iy)=csxy(2,iy);
        csxy(nx,iy)=csxy(nx-1,iy);
    end
    %top and bottom boundary
    for ix = 2:nx-1
        %first derivative of x
        cfx(ix,1)=(C(ix+1,1)-C(ix-1,1))/(2*dx);
        cfx(ix,ny)=(C(ix+1,ny)-C(ix-1,ny))/(2*dx);            
        %first derivative of y
        cfy(ix,1)=(C(ix,2)-C(ix,1))/(dy);
        cfy(ix,ny)=(C(ix,ny)-C(ix,ny-1))/(dy);
        %Second derivative of x,y
        csxy(ix,1)= csxy(ix,2);
        csxy(ix,ny)=csxy(ix,ny-1);
    end
    
    %corners
    
    %first derivative of x
    cfx(1,1)=cfx(1,2); 
    cfx(nx,ny)=cfx(nx,ny-1);
    cfx(nx,1)=cfx(nx,2); 
    cfx(1,ny)=cfx(1,ny-1);
    
    %first derivative of y
    cfy(1,1)=cfy(2,1);
    cfy(nx,1)=cfy(nx-1,1);
    cfy(1,ny)=cfy(2,ny);   
    cfy(nx,ny)=cfy(nx-1,ny);
    
    %second derivative x,y
    csxy(1,1)  =csxy(1,2);
    csxy(nx,1) =csxy(nx,2);
    csxy(nx,ny)=csxy(nx,ny-1);
    csxy(1,ny) =csxy(1,ny-1);
    
    % Reshaping to 1:N
    CFX=reshape(cfx,N,1);
    CFY=reshape(cfy,N,1);
    CSXY=reshape(csxy,N,1);
    C1=reshape(C,N,1);
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%end of concentration derivatives%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    
    %Solution boundary 
    
    for ix = 1:nx
        PSL(ix)=-3*AreaS*C1(ix)/(2*dy);
        PSL1(ix)=2*AreaS*C1(ix)/dy;
        PSL2(ix)=-AreaS*C1(ix)/(2*dy);
    end
     
    % general constants
    
    for ix=1:N
        M1(ix)=C1(ix)/(dy*dy)-CFY(ix)/(2*dy);
        M2(ix)=C1(ix)/(dx*dx)-CFX(ix)/(2*dx);
        M(ix)=-2*C1(ix)/(dx*dx)-2*C1(ix)/(dy*dy);
        M3(ix)=C1(ix)/(dx*dx)+CFX(ix)/(2*dx);
        M4(ix)=C1(ix)/(dy*dy)+CFY(ix)/(2*dy);
        MR(ix)=((UB-UA)/(UB+UA))*CSXY(ix);  
    end
    
    %%%%%%% -nx diagonal
    
    %Variable boundary 1
    Lowdiag2 =[1*ones(nx,1);M1(nx+1:N-nx)';M4(N-nx+1:N-nx+Variable)'+M1(N-nx+1:N-nx+Variable)';zeros((nx-Variable),1);zeros(nx,1)];
    %//////////////////////////

    %%%%%%% -nx+1 diaogonal
    temp0=[0;-1*ones(nx-1,1)];
    Lowdiag1=[temp0;zeros(N,1)];
    
    %%%%%%% -1 diagonal
    temp1=M2'.*repmat([ones(nx-2,1);0;0],[ny,1]);
    temp2=(M2'+M3').*repmat([zeros(nx-2,1);1;0],[ny,1]);
    temp3=temp1+temp2;
    Lowdiag=[zeros(2*nx,1);temp3(nx+1:N)];
    
    %%%%%%% 0 diagonal   
    Diag=[[0;RP(2:nx)'];PSL';M(nx+1:N)'];       
    
    %%%%%%% +1 diagonal
    temp4=M3'.*repmat([0;ones(nx-2,1);0],[ny,1]);
    temp5=(M3'+M2').*repmat([1;zeros(nx-1,1)],[ny,1]);
    temp6=temp4+temp5;
    Updiag=[100;zeros(2*nx,1);temp6(nx+1:N-1)];
    
    %%%%%%% +nx-1 diagonal
    temp7=-1*ones(nx-1,1);
    Updiag1=[zeros(nx,1);temp7;zeros((N-nx+1),1)];
    
    %%%%%%% +nx diagonal
    Updiag2=[zeros(nx,1);1*ones(nx,1);PSL1';M4(nx+1:N-nx)'];
   
    %%%%%%% +2*nx diagonal
    Updiag3=[zeros(3*nx,1);PSL2';zeros(N-3*nx,1)];
    
    Matrix=spdiags([Lowdiag2,Lowdiag1,Lowdiag,Diag,Updiag,Updiag1,Updiag2,Updiag3],[-nx -nx+1 -1 0 1 nx-1 nx 2*nx],N+nx,N+nx); 

    % defining RHS
    RHS0(1)=PPB-Pdel(1);
    for ix =2:nx        
        RHS0(ix)=Pdel(ix-1)-Pdel(ix);  % lower boundary 1
    end
    
    RHS1=AreaS*CFY(1:nx);  
    RHS2=MR(nx+1:N-nx)';
    
    % Variable boundary 1
    RHS3=-MR(N-nx+1:N-nx+(nx-Variable))';
    RHS4=-MR(N-nx+(nx-Variable)+1:N)'-M4(N-nx+(nx-Variable)+1:N)'.*Bpottop((nx-Variable)+1:nx)';    
    RHS=[RHS0';RHS1;RHS2;RHS3;RHS4];
    
    % Solving 
    SOL=Matrix\RHS;

    % Reshaping to nx by ny matrix
    P=reshape(SOL(nx+1:N+nx),nx,ny);  % Solution Potential distribution
    IP=i0*SOL(1:nx);                  % Actual current in the polymer 
    PBoun=SOL(nx+1:2*nx);
    dArea = dx*X0*Width_actual;       % Actual area of polymer between grid points
    
    % Compute actual current density Is
    for ix = 1:nx-1
        IS(ix) = (IP(ix) - IP(ix+1))/dArea; 
    end
    IS(nx) = IS(nx-1);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Part 3 - Finding the Concentration distribution
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % To find the derivative of the potential at each point
    
    for ix = 2:nx-1
        for iy = 2:ny-1
            pfx(ix,iy)=(P(ix+1,iy)-P(ix-1,iy))/(2*dx);
            pfy(ix,iy)=(P(ix,iy+1)-P(ix,iy-1))/(2*dy);
        end
    end
    % left and right boundary
    for iy = 2:ny-1
        % first derivative of ix
        pfx(1,iy)=(P(2,iy)-P(1,iy))/(dx); 
        pfx(nx,iy)=(P(nx,iy)-P(nx-1,iy))/(dx);
        
        % first derivative of iy
        pfy(1,iy)=(P(1,iy+1)-P(1,iy-1))/(2*dy);
        pfy(nx,iy)=(P(nx,iy+1)-P(nx,iy-1))/(2*dy);
        
    end
    % top and bottom boundary
    for ix = 2:nx-1
        %first derivative of x
        pfx(ix,1)=(P(ix+1,1)-P(ix-1,1))/(2*dx);
        pfx(ix,ny)=(P(ix+1,ny)-P(ix-1,ny))/(2*dx);            
        %first derivative of y
        pfy(ix,1)=(P(ix,2)-P(ix,1))/(dy);
        pfy(ix,ny)=(P(ix,ny)-P(ix,ny-1))/(dy);
    end
    % corners
    % first derivative of x
    pfx(1,1)=pfx(1,2); 
    pfx(nx,ny)=pfx(nx,ny-1);
    pfx(nx,1)=pfx(nx,2); 
    pfx(1,ny)=pfx(1,ny-1);
    
    % first derivative of y
    pfy(1,1)=pfy(2,1);
    pfy(nx,1)=pfy(nx-1,1);
    pfy(1,ny)=pfy(2,ny);   
    pfy(nx,ny)=pfy(nx-1,ny);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%% Concentration %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Finding the Alpha and Beta values 
     for ix= 1:nx
        for iy = 1:ny
            k1(ix,iy)=pfx(ix,iy)*charge;
            k2(ix,iy)=pfy(ix,iy)*charge;   
            Gamma1=dt/(2*dx);
            Gamma2=dt/(2*dy);
            
            %Constants for 'X' direction
            if abs(k1(ix,iy))<=1e-12                       
                alpha1(ix,iy)=1/dx;
                beta1(ix,iy) =1/dx;
            else
                alpha1(ix,iy)=-k1(ix,iy)/(1-exp(k1(ix,iy)*dx));
                beta1(ix,iy)=-k1(ix,iy)*exp(k1(ix,iy)*dx)/(1-exp(k1(ix,iy)*dx));                    
            end
    
            %Constants for 'Y' direction
            if abs(k2(ix,iy))<=1e-12                       
                alpha2(ix,iy)=1/dy;
                beta2(ix,iy) =1/dy;
            else
                alpha2(ix,iy)=-k2(ix,iy)/(1-exp(k2(ix,iy)*dy));
                beta2(ix,iy)=-k2(ix,iy)*exp(k2(ix,iy)*dy)/(1-exp(k2(ix,iy)*dy));                  
            end
        end
    end

    %LHS Matrix
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%% n to n+1/2%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % General 
    for iy = 1:ny
        for ix=2:nx-1
            Adiag(ix,iy)=1+Gamma1*(alpha1(ix,iy)+beta1(ix,iy));
            Asubup(ix,iy)=-beta1(ix,iy)*Gamma1;
            Asubdwn(ix,iy)=-alpha1(ix,iy)*Gamma1;
        end
    end
    % boundaries
    for iy=1:ny
        % Left boundary
        Adiag(1,iy)=1+Gamma1*alpha1(ix,1);
        Asubup(1,iy)=-beta1(ix,iy)*Gamma1;
        Asubdwn(1,iy)=0;
        
        % Right boundary  
        Adiag(nx,iy)=1+Gamma1*beta1(ix,ny);
        Asubup(nx,iy)=0;
        Asubdwn(nx,iy)=-alpha1(ix,iy)*Gamma1;
    end

    Diag=reshape(Adiag,1,N)';
    Subdwn=reshape(Asubdwn,1,N)';
    Subup=reshape(Asubup,1,N)';
    
    up=[1000;Subup(1:N-1)];
    dwn=[Subdwn(2:N);1000];
    
    % Matrix
    M=spdiags([dwn,Diag,up],[-1 0 1],N,N);
    
    % RHS n to n+1/2
    
    for iy=2:ny-1
        for ix =1:nx
            R1(ix,iy)=(1-Gamma2*(alpha2(ix,iy)+beta2(ix,iy)))*C(ix,iy)+(beta2(ix,iy+1)*Gamma2)*C(ix,iy+1)+(alpha2(ix,iy-1)*Gamma2)*C(ix,iy-1);
        end
    end

    % top and bottom boundaries
    for ix =1 :nx
        R1(ix,1)= (1-(alpha2(ix,1)*Gamma2))*C(ix,1)+(beta2(ix,2)*Gamma2)*C(ix,2)+(F(ix))*2*Gamma2;
        R1(ix,ny)=(1-(beta2(ix,ny)*Gamma2))*C(ix,ny)+(alpha2(ix,ny-1)*Gamma2)*C(ix,ny-1)-(F(ix))*2*Gamma2;
    end

    % Solving
    R= reshape(R1,1,N);
    
    C1=M\R';
    
    C=reshape(C1,nx,ny);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%% n+1/2 to n+1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % LHS matrix 
    for ix = 1:nx
        for iy=2:ny-1
            Adiag(ix,iy)=1+Gamma2*(alpha2(ix,iy)+beta2(ix,iy));
            Asubup(ix,iy)=-beta2(ix,iy)*Gamma2;
            Asubdwn(ix,iy)=-alpha2(ix,iy)*Gamma2;
        end
    end
    % boundaries
    for ix=1:nx            
        % lower boundary
        Adiag(ix,1)=1+Gamma2*alpha2(ix,1);
        Asubup(ix,1)=0;
        Asubdwn(ix,1)=-alpha2(ix,1)*Gamma2;
        
        % upper boundary
        Adiag(ix,ny)=1+Gamma2*beta2(ix,ny);
        Asubup(ix,ny)=-beta2(ix,ny)*Gamma2;
        Asubdwn(ix,ny)=0;                       
    end

    Diag=reshape(Adiag,1,N)';
    Subdwn=reshape(Asubdwn,1,N)';
    Subup=reshape(Asubup,1,N)';
    
    % Matrix
    M=spdiags([Subdwn,Diag,Subup],[-nx 0 nx],N,N);
    
    % RHS %%%%%%2%%%%%%%%% n+1/2 to%n+1 %%%%%%%%%%%%%%
    
    for iy=1:ny
        for ix =2:nx-1
            R1(ix,iy)=(1-Gamma1*(alpha1(ix,iy)+beta1(ix,iy)))*C(ix,iy)+(beta1(ix+1,iy)*Gamma1)*C(ix+1,iy)+(alpha1(ix-1,iy)*Gamma1)*C(ix-1,iy);                 
        end
    end
    
    for iy =1 :ny
        R1(1,iy)= (1-(alpha1(1,iy)*Gamma1))*C(1,iy)+(beta1(2,iy)*Gamma1)*C(2,iy);            
        R1(nx,iy)=(1-(beta1(nx,iy)*Gamma1))*C(nx,iy)+(alpha1(nx-1,iy)*Gamma1)*C(nx-1,iy);           
    end

    % Solving
    R= reshape(R1,N,1);
    
    C1=M\R;  
    C= reshape(C1,nx,ny);
    
    % Smoothing
    for iy = 4:ny-3
        C(:,iy) = fastsmooth(C(:,iy),5,2,1);
    end

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % End of Part 3
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Part 4 - Finding the current/flux of ions and rate of oxidation
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    for ix=1:nx    
        DPYB(ix) =(4*P(ix,2)-P(ix,3)-3*P(ix,1))/(2*dy);                 % Slope of the potential at the boundary using 3 points (y=1, y=2 and y=3)
        ROX(ix) = IS(ix)/(Mion_actual*Thickness_actual*FC); 
        
        % Limiting ROX when OX -> emax or OX -> emin
        if (OX(ix) + ROX(ix)*dt > emax)
            ROX(ix) = (emax-OX(ix))/dt;
        end
        if (OX(ix) + ROX(ix)*dt < emin)
            ROX(ix) = (emin-OX(ix))/dt;
        end     
        
        % Finding the state of oxidation for the next time step
        OX(ix) = OX(ix)+ROX(ix)*dt;                                       % State of oxidation      
    end

    OX = fastsmooth(OX,5,2,1);

    Cmean = sum(C1)/N;
    OXave = sum(OX)/nx;
    ROXmean=sum(ROX)/nx;
    Pmean=sum(DPYB)/nx;
    
    % To compute polymer currents and potential values
    Ppol(1) = PPB*V0;
    for ix = 2:nx
        Ppol(ix) = Ppol(ix-1) - IP(ix)*RPOLYMER(ix);
    end
    
    % Plot values
    IT(t,1:nx)=IS(1:nx);
    OXPLOT(t,1:nx)=OX(1:nx)/OXmax;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % End of Part 4 
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     
    Pmax(t) = max(IS);      % Storing the maximum current at each time step into an array
    time(t) = dt*T0*t;      % Actual time in seconds
    Actualtime = time(t)
    
    % Display
   figure(1)
   subplot(321)
    mesh(y,x,P*(V0))
    xlabel('y(mm)')
    ylabel('x(mm)')
    title(['Solution Potential (V)    Time = ',num2str(Actualtime),'s'])
    axis([0 YDistance 0 XDistance 0 1.1*Vfinal])
   subplot(322)
    mesh(y,x,C*C0/1000) 
    xlabel('y(mm)')
    ylabel('x(mm)')
    title('Solution Concentration (Moles /liter)')
    %axis([0 YDistance 0 XDistance 0 0.02])
   subplot(323)
    plot(x,IS)
    xlabel('x(mm)')
    ylabel('Current desnity (A/m2)')
    title('Current at the interface')
   axis([0 XDistance 0 70])
   subplot(324)
    plot(x,OX)
    axis([0 XDistance 0 1])
    xlabel('x(mm)')
    ylabel('Degree of oxidation')
    title('Polymer degree of Oxidation ')
  subplot(325)
    plot(x,Ppol)
    axis([0 XDistance 0 1.1*Vfinal])
    xlabel('x(mm)')
    ylabel('Polymer potential (V)')
    title('Polymer potential')
 subplot(326)
    plot(x,IP)
    xlabel('x(mm)')
    ylabel('Polymer current(A)')
    title('Current along the polymer')
  
    IStot = sum(IS)*0.01;
    ISratio = IStot/IP(1);
    IS_int(t) = IStot; 
end

% Max current Vs Time display

   figure(2)
   subplot(221)
    plot(x,OXPLOT(nplot,:),'black',x,OXPLOT(2*nplot,:),'black',x,OXPLOT(3*nplot,:),'black',x,OXPLOT(4*nplot,:),'black',x,OXPLOT(5*nplot,:),'black',x,OXPLOT(6*nplot,:),'black',x,OXPLOT(7*nplot,:),'black',x,OXPLOT(8*nplot,:),'black',x,OXPLOT(9*nplot,:),'black',x,OXPLOT(10*nplot,:),'black')
    xlabel('x(mm)')
    ylabel('Oxidation state(in ratio)')
    title('')
    axis([0 10 0 1])
   subplot(222)
    plot(x,IT(nplot,:),'black',x,IT(2*nplot,:),'black',x,IT(3*nplot,:),'black',x,IT(4*nplot,:),'black',x,IT(5*nplot,:),'black',x,IT(6*nplot,:),'black',x,IT(7*nplot,:),'black',x,IT(8*nplot,:),'black',x,IT(9*nplot,:),'black',x,IT(10*nplot,:),'black')
    xlabel('x(mm)')
    ylabel('Current density(A/m2)')
    title('')
   subplot(223)
     plot(time,IS_int)
     xlabel('t(s)')
     ylabel('Total current (A/m)')
     title('Total current') 
 
end