function f = solve_c(c,params) %returns the residual %unpack the params nz = params.nz;k = params.k; K = params.K; Ha = params.Ha; pa = params.pa; mu = params.mu; rho = params.rho; Da = params.Da; Db = params.Db; Dc = params.Dc; b = params.b ; L = params.L; theta = params.theta; cb_init = params.cb_init;ny = params.ny; y_i = params.y_i; z_i = params.z_i; f = zeros(3*ny*(nz),1); ca = c(1:ny*(nz)); cb = c(ny*(nz)+1:2*ny*(nz)); cc = c(2*ny*(nz)+1:3*ny*(nz)); %convert the array of c's to a matrix ca_matrix = array_to_matrix(ca,ny,nz); cb_matrix = array_to_matrix(cb,ny,nz); cc_matrix = array_to_matrix(cc,ny,nz); %calculate the velocity at each value of y for j=1:ny v(j) = velocity((j-1)*b/ny,params); end %precalculate the reaction rate at every grid point, this saves us the %effort of calculating it every time for every species. for i=1:nz for j=1:ny r(i,j) = reactionRate(ca_matrix(i,j), cb_matrix(i,j), cc_matrix(i,j), params); end end %Equations in the interior of the channel for i = 2:nz-1 for j = 2:ny-1 f((i-1)*ny + j) = Da*doublederivative(ca_matrix,i,j,y_i,2) + Da*doublederivative(ca_matrix,i,j,z_i,1) ... - v(j) * derivative(ca_matrix,i,j,z_i,1) - r(i,j); f((i-1)*ny + j + ny*(nz)) = Db*doublederivative(cb_matrix,i,j,y_i,2) + Db*doublederivative(cb_matrix,i,j,z_i,1) ... - v(j) * derivative(cb_matrix,i,j,z_i,1) - r(i,j); f((i-1)*ny + j + 2*ny*(nz)) = Dc*doublederivative(cc_matrix,i,j,y_i,2) + Dc*doublederivative(cc_matrix,i,j,z_i,1) ... - v(j) * derivative(cc_matrix,i,j,z_i,1) + r(i,j); end end %BC1 z=0 for j=1:ny f(j) = ca_matrix(1,j) - 0; f(j+ny*(nz)) = cb_matrix(1,j) - cb_init; f(j+2*ny*(nz)) = cc_matrix(1,j) - 0; end %BC2 z = L for j=1:ny f((nz-1)*ny+j) = -3*ca_matrix((nz),j) + 4*ca_matrix((nz)-1,j) - ca_matrix((nz)-2,j); f((nz-1)*ny+j + ny*(nz)) = -3*cb_matrix((nz),j) + 4*cb_matrix((nz)-1,j) - cb_matrix((nz)-2,j); f((nz-1)*ny+j +2*ny*(nz)) = -3*cc_matrix((nz),j) + 4*cc_matrix((nz)-1,j) - cc_matrix((nz)-2,j); end %BC y = b for i = 1:nz f((i-1)*ny+ny) = -3*ca_matrix(i,ny) + 4*ca_matrix(i, ny-1) -ca_matrix(i,ny-2); f((i-1)*ny+ny + ny*(nz)) = -3*cb_matrix(i,ny) + 4*cb_matrix(i, ny-1) - cb_matrix(i,ny-2); f((i-1)*ny+ny + 2*ny*(nz)) = -3*cc_matrix(i,ny) + 4*cc_matrix(i, ny-1) - cc_matrix(i,ny-2); end %BC y = 0 for i = 1:nz f((i-1)*ny+1) = ca_matrix(i,1) - pa*Ha; f((i-1)*ny+1+ny*(nz)) = -3*cb_matrix(i,1) + 4*cb_matrix(i,2) - cb_matrix(i,3); f((i-1)*ny+1+2*ny*(nz)) = -3*cc_matrix(i,1) + 4*cc_matrix(i,2) - cc_matrix(i,3); end return;