Problem with fsolve when solving for system of coupled nonlinear equations

I am writing a code for polymerization. The basic equations are diffusion and reaction controlled. They form a set of coupled Partial differential equations. I used finite difference method to convert them into linear and then solve the system using fsolve. Below is the code for the same. It basically consists of the main file and the function file. the function file contains the fsolve syntax. I think I made some mistake there. Please have a look and help me if possible. Perhaps it is the way I have defined my last two equations in the function file

global time_length;
global x_length;

global B0;
global A0;
global A0_bulk;

global const_a;
global const_b;
global const_aa;
global const_aab;
global A0_initial;
global B0_initial;
global const_aaa ;
global K ;
global p;
global hatta;
global delta_t;
global delta_x;
global c0;               % characterstic conc 
global L                 % characterstic length    value is 10^-8
global diffusion_const;  % value is 1.24*10^-11 
global vol_aq;           % value is 110 ml
global vol_org;          % value is 5.5 ml
global reaction_const;   % value is 4*7.92*10^4 m3/kmol s
global area;             % value is 40*10^-4 
global R    ;            % ratio of moles B to A
global V    ;            % ration of volume of org to aq
global nL   ;            % limiting species moles
global dp   ;            % diameter
global m    ;   
global sigmab ;
global sigmac ;
exitflag_vector = zeros(time_length, 1);
%-----------------------------------

time_length = 10;
x_length    = 3;

m = 3;
dp = 1.2*10^-6;
R = 5;             % nb/na
nL = 26.91*10^-4;  % moles
L = 10^-7;
diffusion_const = 1.24*10^-11 ;
vol_aq =  110*10^-6;
vol_org = 5.5*10^-6;
reaction_const = 4*0.85*10^3*10^-3;
area = 6*vol_org/dp ;
V= vol_org/vol_aq;
B0_initial = R*nL/vol_org;
A0_initial = nL/vol_aq;
K = 100;

c0= B0_initial;

hatta = reaction_const*c0*L*L/(diffusion_const)

delta_t = 0.0001;
delta_x = 0.5;

B0  =  zeros(time_length+1, x_length+1);
A0  =  zeros(time_length+1, x_length+1);
A0_bulk = zeros(time_length+1,1);

%%%%non dimension 

A0 = A0/c0;
B0 = B0/c0;
A0_initial = A0_initial/c0;
B0_initial = B0_initial/c0;
A0_bulk = A0_bulk/c0;

const_a = delta_t/(2*(delta_x)^2);
const_b = hatta*delta_t/2;
const_aaa = - area*L/vol_aq;

% initial conditions 

A0_bulk(1) = K*A0_initial;

   A0(1,:)  = 0 ;
   B0(1,1:x_length+1)  = 1 ;
   B0(1, x_length + 1 + 1: (x_length+1)*(m+1)) = 0;               %%should be just m

   % bc for x = 0 and t = 0
   A0(1,1) = A0_bulk(1)/(K);        
   B0(1,1) = B0(1,2)  ;

   for j=1:m                                                      %%should be m-1

   B0(1,(x_length+1)*j+1) = B0(1, (x_length+1)*j+2);

   end

   % for x = infinity and t = 0

   A0(1,x_length+1)  =  0;
   B0(1,x_length+1)  =  B0_initial;    

    for j=2:m+1                                        

   B0(1,(x_length+1)*j) = 0;

   end

    x0 = zeros(((x_length+1)*(m+2))+1,1);               

   for p = 1:time_length

options = optimset('TolFun',1e-6,'MaxFunEvals', 1000000000000,'MaxIter', 1000);

p
[x,feval,exitflag] = fsolve(@myfun_edit2_modified_gud,x0,options);

fprintf('1 Time sample done')
exitflag_vector(p) = exitflag;

A0(p+1,:)    = x(1:x_length+1);
B0(p+1,:)    = x(((x_length+1)*1) +    1 :((x_length+1)*(m+1)) + x_length+1) ;
A0_bulk(p+1) = x(((x_length+1)*(m+2)) + 1);

x0 = x;

end 

%  ....................... function file .........................
  function F = myfun_edit2_modified_gud(x)  
global time_length;
global x_length;

global B0;
global A0;
global MB0;
global MB1;
global MB2;
global A0_bulk;
global B0_bulk;
global sigmab ;
global sigmac ;

global const_a;
global const_b;
global const_aa;
global const_aab;
global A0_initial;
global B0_initial;
global const_aaa ;
global K ;
global p;
global hatta;
global delta_t;
global delta_x;
global c0;               % characterstic conc 
global L                 % characterstic length    value is 10^-8
global diffusion_const;  % value is 1.24*10^-11 
global vol_aq;           % value is 110 ml
global vol_org;          % value is 5.5 ml
global reaction_const;   % value is 4*7.92*10^4 m3/kmol s
global area;             % value is 40*10^-4 
global R   ;              % ratio of moles B to A
global V  ;               % ration of volume of org to aq
global nL;
global m;

m=3;

% Evaluate the vector function
F = zeros((m+2)*(x_length+1)+1,1);                  

A0(p+1,:) =  x(1:x_length+1);
B0(p+1,:) =  x(((x_length+1)*1) +    1 :((x_length+1)*(m+1)) + x_length+1) ;         

A0_bulk(p+1) = x(((x_length+1)*(m+2)) + 1);                                    

j=1:m+1;
%B.C for x = 0
F((x_length-1)*(m+2) + 1)   = A0(p+1,1)   - ((A0_bulk(p+1))/(K)) ;
F((x_length-1)*(m+2) + 1+j) = B0(p+1,(x_length+1)*(j-1)+1) -    B0(p+1,(x_length+1)*(j-1)+2);
F((x_length-1)*(m+2) + m+3) = A0_bulk(p+1) - A0_bulk(p) + ( const_aaa*( A0(p,2)- A0(p,1)));

% Infinity boundary conditions

s=2:m+1;

   F((x_length-1)*(m+2) + m+4)  = A0(p+1,x_length+1)  -   0;
   F((x_length-1)*(m+2) + m+5)  = B0(p+1,x_length+1)  -   B0_initial;
   F((x_length-1)*(m+2) + m+4 + s)  = B0(p+1,(x_length+1)*s) - 0;

for j=1:x_length+1

      for f =1:m+1

      sigmab(p+1,j) = (B0(p+1,(x_length+1)*(f-1)+j)+B0(p+1,(x_length+1)*(f-1)+j));
      sigmab(p,j)  = B0(p,(x_length+1)*(f-1)+j);
      end
  end 

for l =1:m

for j=1:x_length+1

      for f =1:l

      sigmac(p+1,(l-1)*(x_length+1)+j) = (B0(p+1,(x_length+1)*(f-1)+j) + B0(p,(x_length+1)*(f-1)+j))*(B0(p+1,(x_length+1)*(l-f)+j) + B0(p,(x_length+1)*(l-f)+j))/4;
      end
  end
  end

l = 1:m
       q = 2:x_length; 

F(1: (x_length-1) ) = -[B0(p+1,q) - B0(p,q)] + const_a*[B0(p+1,q+1) - 2*B0(p+1,q) + B0(p+1,q-1) + B0(p,q+1) - 2*B0(p,q) + B0(p,q-1)]- [[const_b*[B0(p+1,q)+B0(p,q)]].*[A0(p+1,q)+A0(p,q)]];

F(((x_length-1)*(m+2)+1) :(x_length-1)*(m+3)  ) = -[A0(p+1,q)- A0(p,q)]    + const_a*[A0(p+1,q+1) - 2*A0(p+1,q) + A0(p+1,q-1) + A0(p,q+1) - 2*A0(p,q) + A0(p,q-1)] - [[0.5*const_b*[A0(p+1,q)+ A0(p,q)]].*[sigmab(p+1,q)];

F(((x_length-1) + 1):(x_length-1)*(m+2) ) = -[B0(p+1,(x_length+1)*l+q)- B0(p,(x_length+1)*l+q)] + (const_a/m)*[B0(p+1,(x_length+1)*l+q+1) - 2*B0(p+1,(x_length+1)*l+q) + B0(p+1,(x_length+1)*l+q-1) + B0(p,(x_length+1)*l+q+1) - 2*B0(p,(x_length+1)*l+q) + B0(p,(x_length+1)*l+q-1)] - [[(0.5*const_b)*[A0(p+1,q)+ A0(p,q)]].*[B0(p+1,q)+ B0(p,q)]] + [[(const_b*2)*(A0(p+1,q) + A0(p,q))].*[sigmac(p+1,(x_length+1)*(l-1)+q)]]./[sigmab(p+1,q)] ;   

end

[yanyan li]

i meet the same question.