How to write a functional coefficient f in Parabolic function

I'm solvig a system of 3 nonlinear parabolic PDEs (u1, u2, and u3) on rectangular domain using the parabolic function in matlab. I have some difficulties writing the nonlinear interaction terms. How can I write a function of this form u3/(u2+u1) to be my f1? Any hep is appreciated

xmin=0;xmax=20;ymin=0;ymax=1; 
gd = [3;4;xmin;xmax;xmax;xmin;ymax;ymax;ymin;ymin]; 
sf = 'R1';
ns = [82; 49]; 
[dl, bt] = decsg(gd, sf, ns);
fig1=figure(1);clf
axis([xmin xmax ymin ymax]);
pdegplot(dl, 'edgeLabels', 'on'); 
axis equal
title 'Domain With Edge Labels Displayed'
pg = pdeGeometryFromEdges(dl);
es = pg.Edges([1,2,3,4]);
bc1 = pdeBoundaryConditions(es(1),'q',[0;1;1],'g',[1;20;155]);
bc2 = pdeBoundaryConditions(es(2),'q',[0;0;1],'g',[1;1;48]);
bc3 = pdeBoundaryConditions(es(3),'q',[0;0;1],'g',[1;1;24]);
bc4 = pdeBoundaryConditions(es(4),'q',[0;0;0],'g',[1;1;1]);
[p,e,t] = initmesh(dl,'hmax',0.2); % hmax gives max length of an edge
[p,e,t] = refinemesh(dl,p,e,t);
fig2=figure(2);clf
pdemesh(p,e,t); 
axis equal
title 'Domain With Mesh Displayed'

 np=size(p,2);
 u0s=zeros(np,numberOfPDE);
 u0s(:,1) = ones(np,1); 
 ix = find(and(and(and(p(1,:)>5,p(1,:)<8),p(2,:)>0.2),p(2,:)<=1))
  u0s(ix,1) = zeros(size(ix));

u0s(:,2) = 1*ones(size(np));
u0s(:,3) = 1*ones(size(np));

   u0=zeros(np*numberOfPDE,1);
   u0(1:np)=u0s(:,1);

    u0(np+1:np+np)=u0s(:,1);
    u0(np+1+np:np+np+np)=u0s(:,1);

     tmin=0;tmax=0.1;nframes = 20;
     tlist = linspace(tmin,tmax,nframes);
       numberOfPDE =3;
       pb = pde(numberOfPDE);
       c1 = 1.6;
       c2 = 10.01;
       c3 = 0.1;
       %--------------
        a1 = 0;
        a2 = 0;
        a3 = 0;
        %--------------
        d1 = 0.6;
        d2 = 0.6;
        d3 = 0.1;
        %--------------
        f1 = ??;
        f2 = 0.6;
        f3 = 0.1;

        u1 = parabolic(u0,tlist,pb,p,e,t,c1,a1,f1,d1);
        u2 = parabolic(u0,tlist,pb,p,e,t,c2,a2,f2,d2);
        u3 = parabolic(u0,tlist,pb,p,e,t,c3,a3,f3,d3);

          fig3=figure(3);clf
          colormap(cool);

          for j = 1:nframes,
          pdesurf(p,t,u1(1:np,j));
  axis([xmin xmax ymin ymax 0 1.1]);
  shading interp;
  Mv(j) = getframe;
   end