function [Buff, Cal] = PadeApproximation(N, lambda, mu)
%
%  PadeApproximation(N, lambda, mu) computes the Pade approximation 
% of order N for nondimensionalized Buffer and Ca2+ concentrations for
% given values of parameters (lambda, mu), and returns symbolic 
% functions of nondimensional distance "r" from channel, Buff(r) & Cal(r)
%
%  "Padé approximation of a stationary single-channel Ca2+ nanodomain"
%     Victor Matveev, Sep 10 2016 * Biophysical Journal Supp Info
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

syms An Bn yLow yHigh numLow numHigh denLow denHigh x r eqList coeff Buff Cal;

B = sym('B', [N 1]);  % Symbolic B-coefficient array (denominator coeffs)
A = sym('A', [N 1]);  % Symbolic A-coefficient array (numerator coeffs)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Build the rational functions in variables "x" and "r"="rho"

numHigh = 1;   denHigh = 1;   % Initialize numerator & denominator of x-series
numLow  = r^N; denLow  = r^N; % Initialize numerator & denominator of rho-series

for n = 1 : N
    numHigh = numHigh + A(n)*x^n;    % Build numerator   of x-series
    denHigh = denHigh + B(n)*x^n;    % Build denominator of x-series
    numLow  = numLow + r^(N-n)*A(n); % Build numerator   of rho-series
    denLow  = denLow + r^(N-n)*B(n); % Build denominator of rho-series
end

yHigh = numHigh / denHigh;   % The rational function of "x"
yLow  = numLow  / denLow;    % The rational function of "r"="rho"

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now, plug these rational functions in the main buffer equation, 
% in the rho-form (short-range) and in the x-form (long-range), and Taylor-expand:

eqHigh = taylor( lambda*x^4*diff(yHigh,x,2)                 - x*yHigh  + (1 - yHigh)*(yHigh/mu + 1), 'Order', N+1);
eqLow  = taylor( lambda*(r*diff(yLow,r,2)+2*diff(yLow,r,1)) - yLow + r * (1 - yLow )*( yLow/mu + 1), 'Order', N+1);

eqList = []; n = 1;  % Initialize list of equations for the coefficients

while numel(eqList) < N
    coeff = simplify(subs( diff(eqHigh, x, n), x, 0));  % Taylor coefficient of x-equation
    if coeff ~= 0
        eqList = [eqList, coeff];                       % Add equation term to equation list
    end
    n = n + 1;
end

n = 0;
while numel(eqList) < 2*N
    coeff = simplify(subs( diff(eqLow, r, n), r, 0));  % Taylor coefficient of r-equation
    if coeff ~= 0
        eqList = [eqList, coeff];                      % Add equation term to equation list
    end
    n = n + 1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now, solve the system and pick out solution satisfying all constraints:

S = solve(eqList, 'MaxDegree', 4);  % Find all solutions to the coefficients!
Conds = zeros(numel(S.A1), 4*N);    % Impose 4N conditions on coefficients on each solution:

for n = 1 : N
      An = double( S.(['A', num2str(n)]) );                     % Get all solutions for An
      Bn = double( S.(['B', num2str(n)]) );                     % Get all solutions for Bn
      Conds(:,4*(n-1)+1) = real( An ) > 0;                      % A(n) has to be > 0
      Conds(:,4*(n-1)+2) = real( Bn ) > 0;                      % B(n) has to be > 0
      Conds(:,4*(n-1)+3) = abs( imag(An) ./ real(An) ) < 1e-9;  % A(n) has to be real
      Conds(:,4*(n-1)+4) = abs( imag(Bn) ./ real(Bn) ) < 1e-9;  % B(n) has to be real
end

% Find component of "s" satifying all above physicality conditions:
GoodIndex = find(sum(Conds, 2) == 4*N);  

if (numel(GoodIndex) > 1) 
    fprintf('getPadeApproximation: solution is not unique! \n'); 
    % Not fatal, but contrary to the claim in the paper if true
    GoodIndex = GoodIndex(1);    
end

for n = 1 : N
      An    = real( double( S.(['A', num2str(n)]) ) ); % Plug the satisfactory A(n)
      yHigh = subs(yHigh, A(n), An(GoodIndex) );       % into buffer x-series and
      yLow  = subs(yLow,  A(n), An(GoodIndex) );       % into buffer r-series
      Bn    = real( double( S.(['B', num2str(n)]) ) ); % Plug the satisfactory B(n)
      yHigh = subs(yHigh, B(n), Bn(GoodIndex) );       % into buffer x-series and
      yLow  = subs(yLow,  B(n), Bn(GoodIndex) );       % into buffer r-series 
end

Buff = simplify(yLow);                  % Return the buffer r-series
Cal  = simplify(1/r + (Buff - 1)/mu);   % Use conservation law to find Ca2+

fprintf('>> Pade order N=%d: [Buffer] and [Ca2+] vs distance "r": \n', N);
pretty(Buff)  % Pretty-print the [Buffer] solution
pretty(Cal)   % Pretty-print the [Ca2+]   solution


%%%%%%%%%%%%%%%%%%%%%%%%%%%%  T H E   E N D   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
