function Demo_Convection_Diffusion
    model = createpde;
    L = 0.025; %L
    h = 2E-4; %h
    vb = 1; % Define the correct value.
    gdm = [3;4;0;L;L;0;0;0;h;h];
    g = decsg(gdm,'R1',('R1')');
    geometryFromEdges(model,g);
    generateMesh(model,'Hmax', h/4);
    pdegplot(model,'EdgeLabels','on');            
    D = 2.36E-10;
    specifyCoefficients(model,'m',0,'d',0,'c',D,'a',0,'f',@fcoef);
    applyBoundaryCondition(model,'dirichlet','Edge',4,'u',0);
    applyBoundaryCondition(model,'neumann','Edge',1,'g',@gcoef);
    applyBoundaryCondition(model,'neumann','Edge',3,'g',0);
    R = solvepde(model)
    pdeplot(model,'XYData',R.NodalSolution)
    
    function g = gcoef(location,state)
        Vmax = 0.8E-6;
        Km = 0.457;
        k = Vmax/Km;
        mu = 9.45E-4;
        tauW = 6*mu*vb/h;
        a1 = 9.281E-11;
        a2 = 1.505E-9;
        a3 = 5.624;
        S = a1 + a2*tauW/(a3+tauW);
        g  = k.*state.u - S;
    end
    
    function f = fcoef(location,state)
        yb = location.y/h;
        v = 6*vb*yb.*(1-yb);
        f = v.*state.ux;
    end
end