# MATLAB: Solving a system of PDEs with non constant coefficients and dirichlet and neumann boundary conditions

pdesystem

Hello,
I am trying to solve the following system of equations which deal with the diffusion of chemical species within an electrical double layer of thickness δ:
The diffusion coefficients, with are constant values as well as the reaction rate constants .
The concentration of the species are all dependent on the spacial x coordinate:
with .
For each equation we have one dirichlet condition at (which is a positive value) equal to the initial concentration of species i.
At each equation has a Neumann boundary condition:
I have rewritten the equations in the form:
which is the form that they need to be so I can use the PDE toolbox. I have managed to code the f matrix but I am unclear as to how to write the a matrix, which is a diagonal matrix (see attached images – sorry, couldn't insert them directly here…)
The code runs without issues, but the concentration values do not seem to change along x. I have set the initial conditions for the problem as well, equal to the initial concentrations of the species. I am not sure whether I have set the \$a\$ matrix correctly, but I believe that the f matrix is correct, since I followed the documentation. I tried to do the same with the a matrix, but I was a bit unclear as to how to proceed with it, so I just adapted the example that is given.
Is it possible to solve this kind of system using bvp4c (or bvp5c)? from what I understood of the documentation, these methods apply only to equations with boundary values that do not use derivatives [g(y(a),y(b))=0], right? is there a way to use Neumann boundary conditions with these problems?
Any help will be greatly appreciated!
Thank you all!
I am running the following code:
``%Physical constantsR = 8.3145; %J mol^-1 K^-1T = 298; %KP = 100000; %Pa (1 bar)F = 96485.3329; %s A mol^-1 - Faraday's constantdl = 50*10^-6; %delta [m] - thickness of double layerA = 10^-4; %m^2 - area of electrodeV = 25*10^-6; %m^3 - volume of cellz = 2; %electrons flowingj = [1 5 10 50].*10; %current densities in A m^2%Reaction constantskco2 = 0.84;Kw = 1.00*10^-14; %mol^2 l^-2global k1f, k1f = 2.23*10^3; %l mol^-1 s^-1global k1r, k1r = 1.7*10^-5; %s^-1K1 = k1f/k1r;global k2f, k2f = 6.00*10^9; %l mol^-1 s^-1global k2r, k2r = 7.06*10^4; %l mol^-1 s^-1K2 = k2f/k2r;ka = 1.8*10^-4; %l mol^-1kplus = 1000; %mol m^-3 (1 M) - supporting electrolyte%Diffusion constantsglobal Dco2, Dco2 = 1.67*10^-9; %m^2/sglobal Dhco3, Dhco3 = 1.04*10^-9; %m^2/sglobal Dco32, Dco32 = 8.09*10^-10; %m^2/sglobal Doh, Doh = 5.27*10^-9; %m^2/s%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Determining the initial concentrations of the species%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%CO2 concentration in the gas bubblesco2g = P/(R*T); %mol m^-3co2ini = co2g*kco2; %mol m^-3  -- Initial concentration of CO2%coefficients of the [OH] equation after rearrangingohcoef = [2*K1*K2*co2ini K1*co2ini+1 -kplus -Kw];r = roots(ohcoef);%pick correct rootfor i = 1:length(r)    if (imag(r(i))==0)&&(real(r(i))>0)        ohini = r(i); %mol m^-3  -- Initial concentration of OH    endend  co32ini = K1*K2*ohini^2*co2ini; %mol m^-3 -- Initial concentration of CO32-hco3ini = K1*ohini*co2ini; %mol m^-3 -- Initial concentration of HCO3-inivals = ["[CO2]" "[OH]" "[H]" "[CO32]" "[HCO3]";co2ini ohini hini co32ini hco3ini];%for visualization purpose%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Pseudo steady state concentration profiles in the boundary layer%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%global neqs, neqs = 4;model = createpde(neqs); %4 equations in the systemheight = 5*10^-6; %m = 5 micronsx1 = 0;     y1 = -height/2;x2 = dl;    y2 = -height/2;x3 = dl;    y3 = height/2;x4 = 0;     y4 = height/2;Rect = [3;4;x1;x2;x3;x4;y1;y2;y3;y4];g = decsg(Rect);geom = geometryFromEdges(model,g); %Creting the geometry%{pdegplot(model,'EdgeLabels','on')xlim([0,dl])%}%Neumann boundary conditionsk=1;co2neug = j(k)/(z*F);ohneug = j(k)/(z*F);co32neug = 0;hco3neug = 0;iniconds = [co2ini; hco3ini; co32ini; ohini];eqindex = [1 2 3 4];dirichletconds = [co2ini, hco3ini, co32ini, ohini];neumannconds_g = [co2neug, hco3neug, co32neug, ohneug];%apply Dirichlet boundary condition to edge 2 (right side)applyBoundaryCondition(model,'dirichlet','Edge',2,...                        'u',dirichletconds,'EquationIndex', eqindex);%apply Neumann boundary condition to edge 4 (left side)                    applyBoundaryCondition(model,'neumann','Edge',4,...                        'g',neumannconds_g);ccoef = [Dco2;Dhco3;Dco32;Doh];specifyCoefficients(model, 'm',0,'d',0,'c',ccoef,'a',@acoeffunction,'f', @fcoeffunction);setInitialConditions(model,iniconds);p = generateMesh(model,'Hmax',0.5*10^-6); %0.5micron mesh size%pdeplot(model);results = solvepde(model);u = results.NodalSolution;pdeplot(model,'XYData',u(:,1))%since f is not constant, we must define itfunction f = fcoeffunction(location,state)global k1fglobal k1rglobal k2fglobal k2rglobal neqs; % Number of equationsnr = length(location.x); % Number of columnsf = zeros(neqs,nr); % Allocate f% Now the particular functional form of ff(1,:) = state.u(2,:).*k1r;f(2,:) = state.u(3,:).*k2r + state.u(1,:).*state.u(4,:).*k1f;f(3,:) = state.u(2,:).*state.u(4,:).*k2f;f(4,:) = state.u(2,:).*k1r + state.u(3,:).*k2r;endfunction a = acoeffunction(location,state)global k1fglobal k1rglobal k2fglobal k2rglobal neqs; % Number of equationsnr = length(location.x); % Number of columnsa = zeros(neqs,nr); % Allocate a% Now the particular functional form of aa(1,:) = state.u(4,:).*k1f;a(2,:) = state.u(4,:).*k2f+k1r;a(3,:) = k2r;a(4,:) = state.u(1,:).*k1f + state.u(2,:).*k2f;end``

``function du = f(t,u)CCO2  = u(1);CHCO3 = u(2);CCO3  = u(3);COH   = u(4);du(1:4,1) = u(5:8);du(5,1) = eq1du(5,1) = eq2% ...``
``function res = fb(ua,ub)    % concentration of species at initial momentsres(1) = ub(1) - c1;    res(2) = ub(1) - c2;%...    % derivetivesres(5) = ua(2) - ...%...res = res(:);``