function [Buff Ca] = PadeApproximationQuadraticExplicit(lambda, mu, rr, plotFlag)
%
% PadeApproximationQuadraticExplicit(lambda, mu, rr, plotFlag) computes the 
% 2nd order Pade approximation for nondimensionalized Buffer and Ca2+ 
% concentrations for given values of parameters (lambda, mu), and returns 
% symbolic functions of nondimensional distance "r" from channel, Buff & Ca
%
% If "rr" is an empty list [], it will be populated automatically
% Set plotflag to any non-zero value to see a plot of Buffer concentration
%
%  "Padé approximation of a stationary single-channel Ca2+ nanodomain"
%     Victor Matveev, Sep 10 2016 * Biophysical Journal Supp Info
%
%         Modified and heavily optimized Jan 23, 2019
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

eb = lambda;
q = mu / (1 + mu);

q2  =  q*q;   q3 =  q*q2;  q4  =  q*q3;   
q5  =  q*q4;  q6 =  q*q5;  q7  =  q*q6;
eb2 = eb*eb; eb3 = eb*eb2; eb4 = eb*eb3;

A = 3*(8*eb         + q*(1 - q2));
F =    2*eb*(7 - q) + q*(2 + q - 2*q2 - q3);

D = eb*( 12*eb -  8*q3 + 8*q2        ) + q4*(1 - q)^2;
P = eb*(  8*eb -  4*q3 - 2*q2 + 10*q ) + q2*(1 - q2);
B = eb*( 32*eb - 12*q3 + 2*q2 + 14*q ) + q2*(1 - q2);
E = eb*(  6*eb -    q3 + 3*q2 +  4*q ) + q3*(1 + q - 2*q2);
R = eb*( 24*eb - 10*q3 + 4*q2 +  6*q ) + q3*(1 - q - q2 + q3);

p2 = -243*q6 + 648*q5 - 204*q4 -  594*q3 + 192*q2 +  192*q +  13;
p5 =             6*q5 +   3*q4 -   10*q3 -  10*q2 +    4*q +   7;
p1 =                     10*q4 -    7*q3 -  14*q2 +    7*q +   4;
p3 =                    729*q4 - 1290*q3 -  18*q2 +  569*q +  80;
p6 =                     39*q4 -   69*q3 -  90*q2 +   47*q +  77;
p4 =                                    - 2687*q2 + 2408*q + 832;
p7 =                                       - 9*q2 +    4*q +  13;

G = q7*(q  - 1)^2 +   eb*(2*q5*p1 + eb*q3*p2 +  4*eb2*q2*p3   + 4*eb3*q*p4 + 12000*eb4);
H = q4*(q2 - 1)^2 + 4*eb*(  q3*p5 + eb*q2*p6 + 26*eb2*q *p7 + 400*eb3);

W = B^3 + 6*eb*A*( B*E - 3*eb*A*D + F*(-q*G/3)^(1/2) );

% Final results for the coefficients follow below:

b2 = q * ( W^(1/3) + B + H / W^(1/3) ) / A;
b1 = ( b2*(b2*A - 3*q*P) / (6*eb*q2) - R ) / F;
a1 = b1 - q;
a2 = b2 - q*(b1 - q2);

if numel(rr) < 3
    rr = exp( linspace( log(1e-3), log(100), 80 ) );
end

Buff = real( (rr.^2 + a1*rr + a2) ./ (rr.^2 + b1*rr + b2) );
Ca   = 1./rr + (Buff - 1) ./ mu;

if plotFlag
    plot(rr, Buff);
    title(['2nd Order Pade (\lambda=', num2str(eb), ' \nu=', num2str(1/mu), ')']); %coeffs: ', num2str(real([a1 a2 b1 b2])) ]);
    set(gca, 'xscale', 'log');
end

