Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

MATLAB Chromatography PDE MOL

  1. Oct 4, 2018 #1
    Hello all

    I am using the method of lines to solve the following PDE:

    ## \frac {\partial C} {\partial t} + F\frac {\partial q} {\partial t} + u \frac {dC} {dz} = D_{ax} \frac{\partial^2 C} {\partial z^2} ##

    ## \frac {\partial q} {\partial t} = k (q^{*}-q) ##

    With these initial conditions:
    ## C(0,z) = 0 ##
    ## q(0,z) = 0 ##

    and boundary conditions:
    ## C(t,0) = C_{in} (1) ##
    ## \frac {\partial C} {\partial z}(t,L) = 0 (2) ##

    I discretized the equations above in space and obtained a system of ODE's which I solved with ode15s. To deal with the boundary condition (2) I said:

    ## \frac{C(t,L+1) - C(t,L-1)}{2\Delta z}= 0 ## or ## C(t,L+1) = C(t,L-1) ##

    The problem is that my system gives a solution which tends to explode and makes no sense. My guess is that the way I am discretizing the system is not correct or consistent. I am using between 50-150 nodes.

    My gradient matrix for the 1st derivative (GM1) looks like:

    \begin{bmatrix}
    0 & \frac{1}{2} & 0 & 0 & 0 & 0 & 0 & ... & 0\\
    \frac{-8}{12} & 0 & \frac{8}{12} & \frac{-1}{12} & 0 & 0 & 0 & ... & 0\\
    \frac{1}{12} & \frac{-8}{12} & 0 & \frac{8}{12} & \frac{-1}{12} & 0 & 0 & ... & 0 \\
    0 & \frac{1}{12} & \frac{-8}{12} & 0 & \frac{8}{12} & \frac{-1}{12} & 0 & ... & 0 \\
    ... & ... & ... & ... & ... & ... & ... & ... & ... \\
    0 & ... & 0 & 0 & 0 & \frac{1}{12} & \frac{-8}{12} & \frac{-1}{12} & \frac{8}{12} \\
    0 & ... & 0 & 0 & 0 & 0 & 0 & 0 & 1
    \end{bmatrix}

    and b1 = \begin{bmatrix}
    -\frac{C_{in}}{2} \\
    \frac{C_{in}}{12} \\
    0 \\
    ... \\
    0 \\
    -C_{in}
    \end{bmatrix}

    For the 2nd derivative (GM2):

    \begin{bmatrix}
    -2 & 1 & 0 & 0 & 0 & 0 & 0 & ... & 0\\
    1 & -2 & 1 & 0 & 0 & 0 & 0 & ... & 0\\
    0 & 1 & -2 & 1 & 0 & 0 & 0 & ... & 0 \\
    ... & ... & ... & ... & ... & ... & ... & ... & ... \\
    0 & ... & 0 & 0 & 0 & 0 & 0 & 2 & -2
    \end{bmatrix}

    and b2 = \begin{bmatrix}
    C_{in} \\
    0 \\
    ... \\
    0
    \end{bmatrix}

    In the end I say that:

    ## \frac {\partial C} {\partial z} = \frac{1}{h}GM_{1} \cdot C + \frac{1}{h}b_{1}## and ## \frac {\partial^2 C} {\partial z^2} = \frac{1}{h^2}GM_{2} \cdot C + \frac{1}{h^2}b_{2}##
     
    Last edited: Oct 4, 2018
  2. jcsd
  3. Oct 4, 2018 #2
    My code:

    FUNCTION 1:

    Code (Text):

    function [sol, t, x, C] = LDF_mol(isoType, feedProf, parameter, L, Di, epsb, Q, Cfeed, KLDF, Dax, tpulse, tfinal, opt)
    % Transport-Dispersive Model (TDM) considering mass transfer resistence in the solid to be dominant and
    % using the Linear Driving Force Model (LDF) approach (Glauckauf and Coates, 1947)
    % Change the isotherm model by modifying the isotherm function
    % Change feed profile between pulse (e.g.: chromatografic peak) and step (e.g.: breakthrough experiment)
    % Multicomponent
    % Uses pdepe functin to solve the system of partial differential equations


    global data Grad1 Grad2 res Cs_col

    %% Default arguments (Example)
    if nargin == 0
        isoType =   'mod-langmuir';                                                     % isotherm type, can be 'linear' or 'linear-langmuir'
        feedProf =  'pulse';                                                            % feed profile, can be 'pulse' (e.g.: chromatografic peak) or 'step' (e.g.: breakthrough experiment)
        parameter = [1e3 0.3/1e3 0.05; 1e3 0.1/1e3 0.045];                               % isotherm parameters (depends on the isotherm model chosen)
        L =         10;                                                                 % cm, column length
        Di =        1;                                                                  % cm, column internal diameter
        epsb =      0.4;                                                                % column bulk porosity
        Q =         4;                                                                  % mL/min, flow rate
        Cfeed =    [0.6*300/50 ; 0.6*300/50];                             % g/L, feed concentration
        KLDF =      [60; 60];                                                    % min-1, linear driving force (LDF) mass transfer coefficient
        Dax =       [6.57e-3; 7.57e-3];                                     % cm2/min, axial dispersion coefficient
        tpulse =    0.1;                                                                % min, feed pulse duration. For a step injection set tpulse = tfinal
        tfinal =    5;                                                                  % min, final time for calculation
        opt.npz =   50;                                                                % number of discretization points in z
        opt.npt =   100;                                                                % number of discretization points in t
        opt.fig =   true;                                                               % true - show figures; false - do not show figures

        Cs_wash = 0; % Salt concentration in mM
        twash = 1;
     
        t0_elu = tpulse+twash; telu = 8; Cf = 80;
        Cs_elu = @(t) 0+(Cf-0)/telu*(t-t0_elu);

        Cs_reg = 80;
        treg = 3;
     
        Cs_eq = 0;
        teq = 3;
    end

    %% Calculations
    addpath('../')
    data = struct('isoType',isoType,'feedProf',feedProf,'parameter',parameter,'Cfeed',Cfeed,'Cs_wash',Cs_wash,'Cs_elu',Cs_elu,'Cs_reg',Cs_reg,'Cs_eq',Cs_eq,'npz',opt.npz,'npt',opt.npt,'Q',Q,'L',L,'epsb',epsb,'tpulse',tpulse,'twash',twash,'t0_elu',t0_elu,'telu',telu,'treg',treg,'teq',teq,'KLDF',KLDF,'Dax',Dax);
    nc = length(Cfeed);
    data.nc = nc;
    A = pi()*Di^2/4; % cm2
    data.F = (1-epsb)/epsb;
    data.ui = Q/A; % cm/min

    % Run MOL
    t = linspace(0,tfinal,opt.npt);
    z = linspace(0,L,opt.npz);

    [Grad1,Grad2] = GradientMatrix(opt.npz);
    u0 = zeros(2*nc*opt.npz,1);

    [t,uout] = ode15s(@odefun, t, u0);
    for CompNo = 1:data.nc
        idx1 = 1+(CompNo-1)*opt.npz; idx2 = CompNo*opt.npz;
     
        C{CompNo} = uout(:,idx1:idx2);
    end
    res.t = t'; res.z = z'; res.C = C;

    Cs_col = zeros(opt.npt,opt.npz);

    Cs_col(:,1) = Cs_profile(t);
    for i = 1:eek:pt.npt
        for j = 2:eek:pt.npz
            if t(i)-z(j)/data.ui >= 0
                Cs_col(i,j) = Cs_profile(t(i)-z(j)/data.ui);
            else
                Cs_col(i,j) = 0;
            end
        end
    end

    ymat = [];
    for CompNo = 1:data.nc
        ymat = [ymat C{CompNo}(:,end)];
    end  
    [ax,h1,h2] = plotyy(t,ymat,t,Cs_col(:,1));
    axis([0 tfinal  0 inf]) % fix the axes
    set(ax(1),'Box','off'); set(ax(1),'YTick',0:0.2:2); set(ax(2),'YTick',0:20:100,'YLim',[0 100]);
    xlabel('t')
    ylabel('C')

    % Concentration over time with slider figure
        figure
        hold all
        slidermin = 1; % t = 0
        slidermax = opt.npt; % t = tfinal
        ymax_aux = [];
        for CompNo = 1:data.nc-1
            res.hslider(CompNo) = plot(z',C{CompNo}(end,:));
            ymax_aux = [ymax_aux max(C{CompNo}(2,:))];
        end
        ymax_aux = [ymax_aux max(C{data.nc}(2,:))];
        res.ymax = str2double(num2str(max(ymax_aux),1));
        axis([0 L 0 2]);
     
        [ax,h1,h2] = plotyy(z',C{data.nc}(end,:),z',Cs_col(end,:));
        set(ax(1),'Box','off'); set(ax(2),'YTick',0:20:100,'YLim',[0 100]);
        hold off
        res.hslider(data.nc) = h1;
        res.hslider(data.nc+1) = h2;
     
        xlabel('Position, z')
        ylabel('Concentration, C')
        anno_text = sprintf('t = %.1f', 0);
        res.time_box = annotation('textbox',[.55 .48 .4 .5],'String',anno_text,'EdgeColor', 'none','FitBoxToText','on');
        sld = uicontrol('Style', 'slider',...
            'Min',slidermin,'Max',slidermax,'Value',slidermin,...
            'SliderStep',[1 1]./(slidermax-slidermin),...
            'Position', [390 390 120 20],...
            'Callback', @sliderplot);

    %% Auxiliary functions
    function sliderplot(source,event)
    % Handles the UI slider in the Concentration vs Position plot
    global data res Cs_col

    val = round(get(source,'Value'));
    for CompNo = 1:data.nc
        set(res.hslider(CompNo),'YData',res.C{CompNo}(val,:));
    end
    set(res.hslider(data.nc+1),'YData',Cs_col(val,:));

    anno_text = sprintf('t = %.2f',res.t(val));
    set(res.time_box,'String',anno_text);

    function udot = odefun(t,u)
        global data Grad1 Grad2
     
        nc = data.nc; npz = data.npz; L = data.L;
        z = linspace(0,L,npz); h = L/(npz-1);
     
        Cin = setFeedProfile(data.feedProf, t, data.tpulse, data.Cfeed);
     
        C = u(1:nc*npz);
        for CompNo = 1:nc
            idx1 = 1+(CompNo-1)*npz;
            C(idx1) = Cin(CompNo);
        end
     
        q = u(nc*npz+1:end);
        Cs = zeros(npz,1);
     
        udot = zeros(2*nc*npz,1);  
        for j = 1:npz
            if t-z(j)/data.ui >= 0
                Cs(j) = Cs_profile(t-z(j)/data.ui);
            else
                Cs(j) = 0;
            end
        end
     
        qeq = isotherm1(data.isoType, C, Cs, nc, data.parameter);

        dCdz = zeros(nc*npz,1); d2Cdz2 = zeros(nc*npz,1);
        for CompNo = 1:nc
            idx1 = 1+(CompNo-1)*npz; idx2 = CompNo*npz;
            idx3 = npz*(nc+CompNo-1)+1; idx4 = (nc+CompNo)*npz;
         
            b1 = zeros(npz,1); b1(2) = -Cin(CompNo)/2; b1(3) = Cin(CompNo)/12;
            b1(end) = -C(idx2);
            b2 = zeros(npz,1); b2(2) = Cin(CompNo);

            dCdz(idx1:idx2) = Grad1*C(idx1:idx2)./h + b1./h;
            d2Cdz2(idx1:idx2) = Grad2*C(idx1:idx2)./h^2 + b2./h^2;
         
            udot(idx3:idx4) = data.KLDF(CompNo)*(qeq(:,CompNo)-q(idx1:idx2));
            udot(idx1:idx2) = -data.ui*dCdz(idx1:idx2) - data.F*udot(idx3:idx4) + data.Dax(CompNo)*d2Cdz2(idx1:idx2);

            udot(idx1)=0; udot(idx3) = 0;
        end
     
    %% Feed profile
    function Cinj=setFeedProfile(feedProf, t, tpulse, Cfeed)
    % Defines the feed profile. Can be 'pulse' (e.g.: chromatografic peak) or 'step' (e.g.: breakthrough experiment)

    if strcmp(feedProf,'pulse')
        if t<=tpulse
            Cinj = Cfeed;
        else
            Cinj = zeros(1,length(Cfeed));
        end
     
    elseif strcmp(feedProf,'step')
        Cinj = Cfeed;
     
    else
        error('Invalid feed profile. feedProf must be "step" or "pulse"')
    end

    function Cs_in = Cs_profile(t)
    global data

    Cs_in = zeros(length(t),1);

    for i = 1:length(t)
     
        if t(i) <= data.tpulse
            Cs_in(i) = 0;
        elseif t(i) <= data.tpulse+data.twash
            Cs_in(i) = data.Cs_wash;
        elseif t(i) <= data.tpulse+data.twash+data.telu
            Cs_in(i) = data.Cs_elu(t(i));
        elseif t(i) <= data.tpulse+data.twash+data.telu+data.treg
            Cs_in(i) = data.Cs_reg;
        else
            Cs_in(i) = data.Cs_eq;
        end
    end

    function [Grad1,Grad2] = GradientMatrix(N)

    Grad1 = zeros(N,N);
    Grad1(2,2) = 1/2; Grad1(3,1:4) = [-8/12, 0, 8/12, -1/12];
    Grad1(end-1,end-3:end) = [1/12 -8/12 -1/12 8/12]; Grad1(end,end) = 1;

    for i = 4:N-2
        Grad1(i,i-3:i+1) = [1/12, -8/12, 0, 8/12, -1/12];
    end

    Grad2 = zeros(N,N); Grad2(2,1:2) = [-2, 1];
    Grad2(end,end-1:end) = [2,-2];

    for i = 3:N-1
        Grad2(i,i-2:i) = [1, -2, 1];
    end
     
     
  4. Oct 4, 2018 #3
    FUNCTION 2

    Code (Text):

    %% Isotherm Library

    function q=isotherm1(isoType, C, Cs, nc, parameter)
    % Defines the isotherm
    parameter = transpose(parameter); N = length(Cs);

    q = zeros(N,nc);
    Cmat = zeros(N,nc);
    for i = 1:nc
        idx1 = 1+(i-1)*N; idx2 = i*N;
        Cmat(:,i) = C(idx1:idx2);
    end

    if strcmp(isoType,'mod-langmuir')
    % Linear-Langmuir isotherm:   q1(C1,C2) = a1*b1*exp(-m1*Cs)*C1/(1+b1*exp(-m1*Cs)*C1) , parameter = [a1 ; b1 ; m1]
        for i = 1:N
            langmuirDenominator = 1;
            for j = 1:nc
                langmuirDenominator = langmuirDenominator +  parameter(2,j)*exp(-parameter(3,j)*Cs(i))*Cmat(i,j);
            end
            q(i,:) = parameter(1,:).*parameter(2,:).*exp(-parameter(3,:)*Cs(i)).*Cmat(i,:)/langmuirDenominator;
        end
     
    else
        error('Invalid isotherm type. isoType must be "linear" or "langmuir" or "linear-langmuir"')
    end
     
     
  5. Oct 4, 2018 #4
    I think now I kind of solved it. However when the wave reaches the right boundary there is some noise in the wave. Do you know what could be the reason?

    Code (Text):

    function [sol, t, x, C] = LDF_mol(isoType, feedProf, parameter, L, Di, epsb, Q, Cfeed, KLDF, Dax, tpulse, tfinal, opt)
    % Transport-Dispersive Model (TDM) considering mass transfer resistence in the solid to be dominant and
    % using the Linear Driving Force Model (LDF) approach (Glauckauf and Coates, 1947)
    % Change the isotherm model by modifying the isotherm function
    % Change feed profile between pulse (e.g.: chromatografic peak) and step (e.g.: breakthrough experiment)
    % Multicomponent
    % Uses pdepe functin to solve the system of partial differential equations


    global data Grad1 Grad2 res Cs_col

    %% Default arguments (Example)
    if nargin == 0
        isoType =   'mod-langmuir';                                                     % isotherm type, can be 'linear' or 'linear-langmuir'
        feedProf =  'pulse';                                                            % feed profile, can be 'pulse' (e.g.: chromatografic peak) or 'step' (e.g.: breakthrough experiment)
        parameter = [1e3 1.5/1e3 0.15; 1e3 2/1e3 0.045];                               % isotherm parameters (depends on the isotherm model chosen)
        L =         10;                                                                 % cm, column length
        Di =        1;                                                                  % cm, column internal diameter
        epsb =      0.4;                                                                % column bulk porosity
        Q =         4;                                                                  % mL/min, flow rate
        Cfeed =    [0.6*300/50 ; 0.6*300/50];                             % g/L, feed concentration
        KLDF =      [60; 60];                                                    % min-1, linear driving force (LDF) mass transfer coefficient
        Dax =       [6.57e-3; 7.57e-3];                                     % cm2/min, axial dispersion coefficient
        tpulse =    0.1;                                                                % min, feed pulse duration. For a step injection set tpulse = tfinal
        tfinal =    8;                                                                  % min, final time for calculation
        opt.npz =   300;                                                                % number of discretization points in z
        opt.npt =   100;                                                                % number of discretization points in t
        opt.fig =   true;                                                               % true - show figures; false - do not show figures

        Cs_wash = 0; % Salt concentration in mM
        twash = 1;
       
        t0_elu = tpulse+twash; telu = 8; Cf = 80;
        Cs_elu = @(t) 0+(Cf-0)/telu*(t-t0_elu);

        Cs_reg = 80;
        treg = 3;
       
        Cs_eq = 0;
        teq = 3;
    end

    %% Calculations
    addpath('../')
    data = struct('isoType',isoType,'feedProf',feedProf,'parameter',parameter,'Cfeed',Cfeed,'Cs_wash',Cs_wash,'Cs_elu',Cs_elu,'Cs_reg',Cs_reg,'Cs_eq',Cs_eq,'npz',opt.npz,'npt',opt.npt,'Q',Q,'L',L,'epsb',epsb,'tpulse',tpulse,'twash',twash,'t0_elu',t0_elu,'telu',telu,'treg',treg,'teq',teq,'KLDF',KLDF,'Dax',Dax);
    nc = length(Cfeed);
    data.nc = nc;
    A = pi()*Di^2/4; % cm2
    data.F = (1-epsb)/epsb;
    data.ui = Q/A; % cm/min

    % Run MOL
    t = linspace(0,tfinal,opt.npt);
    z = linspace(0,L,opt.npz);

    [Grad1,Grad2] = GradientMatrix(opt.npz);
    u0 = zeros(2*nc*opt.npz,1);

    [t,uout] = ode15s(@odefun, t, u0);
    for CompNo = 1:data.nc
        idx1 = 1+(CompNo-1)*opt.npz; idx2 = CompNo*opt.npz;
       
        C{CompNo} = uout(:,idx1:idx2);
    end
    res.t = t'; res.z = z'; res.C = C;

    Cs_col = zeros(opt.npt,opt.npz);

    Cs_col(:,1) = Cs_profile(t);
    for i = 1:opt.npt
        for j = 2:opt.npz
            if t(i)-z(j)/data.ui >= 0
                Cs_col(i,j) = Cs_profile(t(i)-z(j)/data.ui);
            else
                Cs_col(i,j) = 0;
            end
        end
    end

    ymat = [];
    for CompNo = 1:data.nc
        ymat = [ymat C{CompNo}(:,end)];
    end  
    [ax,h1,h2] = plotyy(t,ymat,t,Cs_col(:,1));
    axis([0 tfinal  0 inf]) % fix the axes
    set(ax(1),'Box','off'); set(ax(1),'YTick',0:0.2:2); set(ax(2),'YTick',0:20:100,'YLim',[0 100]);
    xlabel('t')
    ylabel('C')

    % Concentration over time with slider figure
        figure
        hold all
        slidermin = 1; % t = 0
        slidermax = opt.npt; % t = tfinal
        ymax_aux = [];
        for CompNo = 1:data.nc-1
            res.hslider(CompNo) = plot(z',C{CompNo}(end,:));
            ymax_aux = [ymax_aux max(C{CompNo}(2,:))];
        end
        ymax_aux = [ymax_aux max(C{data.nc}(2,:))];
        res.ymax = str2double(num2str(max(ymax_aux),1));
        axis([0 L 0 2]);
       
        [ax,h1,h2] = plotyy(z',C{data.nc}(end,:),z',Cs_col(end,:));
        set(ax(1),'Box','off'); set(ax(2),'YTick',0:20:100,'YLim',[0 100]);
        hold off
        res.hslider(data.nc) = h1;
        res.hslider(data.nc+1) = h2;
       
        xlabel('Position, z')
        ylabel('Concentration, C')
        anno_text = sprintf('t = %.1f', 0);
        res.time_box = annotation('textbox',[.55 .48 .4 .5],'String',anno_text,'EdgeColor', 'none','FitBoxToText','on');
        sld = uicontrol('Style', 'slider',...
            'Min',slidermin,'Max',slidermax,'Value',slidermin,...
            'SliderStep',[1 1]./(slidermax-slidermin),...
            'Position', [390 390 120 20],...
            'Callback', @sliderplot);

    %% Auxiliary functions
    function sliderplot(source,event)
    % Handles the UI slider in the Concentration vs Position plot
    global data res Cs_col

    val = round(get(source,'Value'));
    for CompNo = 1:data.nc
        set(res.hslider(CompNo),'YData',res.C{CompNo}(val,:));
    end
    set(res.hslider(data.nc+1),'YData',Cs_col(val,:));

    anno_text = sprintf('t = %.2f',res.t(val));
    set(res.time_box,'String',anno_text);

    function udot = odefun(t,u)
        global data Grad1 Grad2
       
        nc = data.nc; npz = data.npz; L = data.L;
        z = linspace(0,L,npz); h = L/(npz-1);
       
        Cin = setFeedProfile(data.feedProf, t, data.tpulse, data.Cfeed);
       
        C = u(1:nc*npz);
        for CompNo = 1:nc
            idx1 = 1+(CompNo-1)*npz;
            C(idx1) = Cin(CompNo);
        end
       
        q = u(nc*npz+1:end);
        Cs = zeros(npz,1);
       
        udot = zeros(2*nc*npz,1);  
        for j = 1:npz
            if t-z(j)/data.ui >= 0
                Cs(j) = Cs_profile(t-z(j)/data.ui);
            else
                Cs(j) = 0;
            end
        end
       
        qeq = isotherm1(data.isoType, C, Cs, nc, data.parameter);

        dCdz = zeros(nc*npz,1); d2Cdz2 = zeros(nc*npz,1);
        for CompNo = 1:nc
            idx1 = 1+(CompNo-1)*npz; idx2 = CompNo*npz;
            idx3 = npz*(nc+CompNo-1)+1; idx4 = (nc+CompNo)*npz;
           
            b1 = zeros(npz-1,1); b1(1) = -Cin(CompNo)/2; b1(2) = Cin(CompNo)/12;
            b1(end) = -C(idx2);
            b2 = zeros(npz-1,1); b2(1) = Cin(CompNo);
           
            dCdz(idx1) = 3.14; % Value does not matter
            d2Cdz2(idx1) = 3.14; % Value does not matter
           
            dCdz(idx1+1:idx2) = Grad1*C(idx1+1:idx2)./h + b1./h;
            d2Cdz2(idx1+1:idx2) = Grad2*C(idx1+1:idx2)./h^2 + b2./h^2;
           
            udot(idx3:idx4) = data.KLDF(CompNo)*(qeq(:,CompNo)-q(idx1:idx2));
            udot(idx1:idx2) = -data.ui*dCdz(idx1:idx2) - data.F*udot(idx3:idx4) + data.Dax(CompNo)*d2Cdz2(idx1:idx2);

            udot(idx1)=0; udot(idx3) = 0;
        end
       
    %% Feed profile
    function Cinj=setFeedProfile(feedProf, t, tpulse, Cfeed)
    % Defines the feed profile. Can be 'pulse' (e.g.: chromatografic peak) or 'step' (e.g.: breakthrough experiment)

    if strcmp(feedProf,'pulse')
        if t<=tpulse
            Cinj = Cfeed;
        else
            Cinj = zeros(1,length(Cfeed));
        end
       
    elseif strcmp(feedProf,'step')
        Cinj = Cfeed;
       
    else
        error('Invalid feed profile. feedProf must be "step" or "pulse"')
    end

    function Cs_in = Cs_profile(t)
    global data

    Cs_in = zeros(length(t),1);

    for i = 1:length(t)
       
        if t(i) <= data.tpulse
            Cs_in(i) = 0;
        elseif t(i) <= data.tpulse+data.twash
            Cs_in(i) = data.Cs_wash;
        elseif t(i) <= data.tpulse+data.twash+data.telu
            Cs_in(i) = data.Cs_elu(t(i));
        elseif t(i) <= data.tpulse+data.twash+data.telu+data.treg
            Cs_in(i) = data.Cs_reg;
        else
            Cs_in(i) = data.Cs_eq;
        end
    end

    function [Grad1,Grad2] = GradientMatrix(N)

    % Grad1 = zeros(N,N);
    % Grad1(2,2) = 1/2; Grad1(3,1:4) = [-8/12, 0, 8/12, -1/12];
    % Grad1(end-1,end-3:end) = [1/12 -8/12 -1/12 8/12]; Grad1(end,end) = 1;
    %
    % for i = 4:N-2
    %     Grad1(i,i-3:i+1) = [1/12, -8/12, 0, 8/12, -1/12];
    % end
    Grad1 = zeros(N-1,N-1);
    Grad1(1,2) = 1/2; Grad1(2,1:4) = [-8/12, 0, 8/12, -1/12];
    Grad1(end-1,end-3:end) = [1/12 -8/12 -1/12 8/12]; Grad1(end,end) = 1;

    for i = 3:N-3
        Grad1(i,i-2:i+2) = [1/12, -8/12, 0, 8/12, -1/12];
    end

    Grad2 = zeros(N-1,N-1); Grad2(1,1:2) = [-2, 1];
    Grad2(end,end-1:end) = [2,-2];

    for i = 2:N-2
        Grad2(i,i-1:i+1) = [1, -2, 1];
    end
     
     
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Have something to add?