How can make this code converge to have the values

function dcdt=anode(tfn,cfn);

global j x ccc in phi1 phi2 etas1 v phi10;

%===============CALCULATE IN===========================

ccc=cfn';

v=0.0000001;

current=0.001;   % Current in A/m2

d=1e-9;          % m2/s

jflux=current;   % Flux for Fick's Law

epsilon=0.5;     % Porosity

tn=0.33;         % Transfer Number

f=96487;         % Coulombs/equivalent

ka1=2e-6;        % m/s

alphac1=0.5;

kc1=2e-6;        % m/s

alphaa1=0.5;
cmax=1000;       % mol/m3

u1=0;           
%r=83.14;          %m3 * Pa / mol / K

r=8.314;

t=298.14;        %K

k1eff=32.3;       %S/m

k2eff=0.06;       %S/m 

a=188400;            % Interfacial Area per Volume

lcfn=log(cfn);

for i=3:(j-2) ;       %define second derivatives, do not use border values.

      ddlcidx(i)=(lcfn(i-1)-2*lcfn(i)+lcfn(i+1))/(x(i+1)-x(i))^2;
  end
  ddlcidx(2)=((lcfn(3)-lcfn(2))/(x(i+1)-x(i))+0)/(x(i+1)-x(i));

ddlcidx(j-1)=((lcfn(j-2)-lcfn(j-1))/(x(i+1)-x(i))+0)/(x(i+1)-x(i));

ddlcidx(j)=((lcfn(j-1)-lcfn(j))/(x(i+1)-x(i))+0)/(x(i+1)-x(i));

for i=1:j;

     etas1(i)=phi1(i)-phi2(i)-u1;

     i0(i)=f*ka1^alphac1*kc1^alphaa1*(cmax-cfn(i))^alphac1*cfn(i)^alphaa1;

     in(i)=i0(i)*(exp(alphaa1*f*etas1(i)/r/t)-exp(-alphac1*f*etas1(i)/r/t));
  end;

in=real(in);

%=========CALCULATE IN AND VOLTAGES=======

phi1(1)=0;               % BC Solid Voltage set solid as ground.

phi2(1)=0;               % BC    ???????  to get program to run.

dphi1(1)=-jflux/k1eff;    % BC Solid Voltage, specified current at collector.

dphi2(1)=-0.01;              % BC Electrolyte Voltage, no current at wall. 

ddphi1(1)=a*in(1)/k1eff;

ddphi2(1)=(2*r*t*k2eff/f*0.67*ddlcidx(1)-a*in(1))/k2eff;

for i=2:j;

     ddphi1(i)=a*in(i)/k1eff;

     ddphi2(i)=(2*r*t*k2eff/f*0.67*ddlcidx(i)-a*in(i))/k2eff;

     dphi1(i)=ddphi1(i)*(x(i)-x(i-1))+dphi1(i-1);

     dphi2(i)=ddphi2(i)*(x(i)-x(i-1))+dphi2(i-1);

     phi1(i)=dphi1(i)*(x(i)-x(i-1))+phi1(i-1);

     phi2(i)=dphi2(i)*(x(i)-x(i-1))+phi2(i-1);
  end;
  %=========NOW CALCULATE CONCENTRATIONS==================================

dcidx(1)=0;           %BC of Fick's Law and constant current.  OPTIONAL.

dcidx(j)=-jflux/d;    %BC of Fick's Law and constant current.  OPTIONAL.

for i=3:(j-2) ;       %define second derivatives, do not use border values.

      ddcidx(i)=(cfn(i-1)-2*cfn(i)+cfn(i+1))/(x(i+1)-x(i))^2;
  end
  ddcidx(2)=((cfn(3)-cfn(2))/(x(i+1)-x(i))+dcidx(1))/(x(i+1)-x(i));

ddcidx(j-1)=((cfn(j-2)-cfn(j-1))/(x(i+1)-x(i))+dcidx(j))/(x(i+1)-x(i));

ddcidx(j)=ddcidx(j-1); %OPTIONAL

%dcdt=d/epsilon*ddcidx+q(1+tn)*a/f*in;

dcdt=d/epsilon*ddcidx-v*dcidx+(1+tn)*a/f*in;

dcdt(1)=dcdt(2);       %OPTIONAL

dcdt(j)=dcdt(j-1);     %OPTIONAL

dcdt=dcdt';

%===============SET UP VECTOR VARIABLES================

global j x ccc in phi1 phi2 etas1 phi10 cc;

%=================INITIAL CONDITIONS===================

j=101;              % x=0 corresponds to anode (Li --> Li+ + e) wall.

                  % x=1 corresponds to separator.

c0=999;

x(1)=0;

c(1)=c0;

for i=2:j;

      c(i)=c0;

      x(i)=x(i-1)+0.01/(j-1);   % 0.01 meters

end;

phi1=zeros(1,j);

phi1=real(phi1);

phi2=zeros(1,j);

phi2=real(phi2);

cc=real(cc);

tspan=[0 100];

%========CALL ON ODE45 TO CALCULATE CONCENTRATIONS AND VOLTAGES=========

[tt,cc]=ode45('anode',tspan,c);

[Walter Roberson]

please reformat so that people can read this.

http://www.mathworks.com/matlabcentral/answers/13205-tutorial-how-to-format-your-question-with-markup

[Jan]

There is a bunch of global variables. Even for the author itself this is hard to debug.

Do you have any idea, why this program does not do what you expect? Could you describe, what you expect and what you see instead? Do you assume a bug in the code or is the model not sufficient?

What should we assume when reading this line:

phi2(1)=0; % BC ??????? to get program to run.