a1 = (-Vw*dt/(2*dz) + D*dt/(dz*dz))/Rf;

a2 = 1 - 2*D*dt/(Rf*dz*dz);

a3 = (Vw*dt/(2*dz) + D*dt/(dz*dz))/Rf;

c1 = -Vw*dt/(2*dz) - V*dt/(2*dz) + tau*D*dt/(dz*dz);

c2 = 1 - 2*tau*D*dt/(dz*dz);

c3 = Vw*dt/(2*dz) + V*dt/(2*dz) + tau*D*dt/(dz*dz);

b1 = mum/((ks + S(i,j-1))*Y);

b2 = (M(i,j-1) + rhos*Ma(i,j-1)/theta)*dt/Rf;

S(i,j) = S(i+1,j-1)*a1 + S(i,j-1)*(a2-b1*b2) + S(i-1,j-1)*a3;

d1 = ka*(1 - Ma(i,j-1)/Mam)*dt;

d2 = (mum*S(i,j-1)/(ks + S(i,j-1)) - b)*dt;

d3 = (ky*rhos*dt/theta)*Ma(i,j-1);

M(i,j) = M(i+1,j-1)*c1 + M(i,j-1)*(c2-d1+d2) + M(i-1,j-1)*c3 + d3;

b1 = mum/((ks + S(i,j-1))*Y);

b2 = (M(i,j-1) + rhos*Ma(i,j-1)/theta)*dt/Rf;

S(i,j) = S(i-1,j-1)*a1 + S(i,j-1)*(a2-b1*b2) + S(i-1,j-1)*a3;

d1 = ka*(1 - Ma(i,j-1)/Mam)*dt;

d2 = (mum*S(i,j-1)/(ks + S(i,j-1)) - b)*dt;

d3 = ky*rhos*dt*Ma(i,j-1);

M(i,j) = M(i-1,j-1)*c1 + M(i,j-1)*(c2-d1+d2) + M(i-1,j-1)*c3 + d3;

e1 = 1 - ka*theta*M(i,j-1)/(rhos*Mam);

e2 = mum*S(i,j-1)/(ks + S(i,j-1)) - b;

e3 = ka*theta*dt*M(i,j-1)/rhos;

Ma(i,j) = Ma(i,j-1)*(e1 - ky + e2)*dt + e3;