I don't have enough reputation to write a comment, so this one goes down as an answer.
I take it you've written a code to solve the angular distribution $\psi$ numerically for a given concentration $c$, using the self-consistent procedure outlined in the paper. Now that you can write the free energy as a function of phi, and the pressure is the $c$-derivative of the free energy, you can, in principle, get the pressure for any value of $c$.
I am slightly confused by the fact that it seems, in the derivation for pressure, as if the $\psi$s do not depend on $c$. And they obviously do. So I believe there is an implicit approximation here. If you keep to the same approximation, you can write down the derivative of the pressure by hand. You have then a way to evaluate both the pressure and its derivative at any given $c$ (this is why I am not sure what you mean when you say that you only know them for certain values of $c$). Now that you have these functions (and similarly for the chemical potentials), writing down the Newton-Raphson scheme should pose no problem (you have a two variable, two function setting, so you'll have to differentiate both functions with respect to each variable, and then invert the resulting Jacobian matrix).
So I decided to write the code myself, in "Python". I say "Python" instead of Python, for I used it as a quick drop in toolbox. I uses ipython with the -pylab switch so it'll load parts of scipy, numpy for computations.
Ntheta = 128
ks = arange(Ntheta+1)[1:]
thetas = pi/2.*ks/(Ntheta+1)
deltas = zeros_like(thetas)
deltas = 1 - (cos(thetas) + cos(thetas))/2.
deltas[1:-1] = (cos(thetas[:-2]) - cos(thetas[2:]))/2.
deltas[-1] = (cos(thetas[-1]) + cos(thetas[-2]))/2.
Nphi = 256
js = arange(Nphi+1)[1:]
phis = 2.*pi*js/(Nphi+1)
Ks = zeros((Ntheta, Ntheta))
for k in range(Ntheta):
for l in range(Ntheta):
g = sqrt(1. - (cos(thetas[k])*cos(thetas[l]) + sin(thetas[k])*sin(thetas[l])*cos(phis))**2)
Ks[k, l] = 2.*pi/(Nphi+1) * (3./2*g + sum(g[1:-1]) + 3./2*g[-1])
taus = (c/pi)**2 * exp(-2*c**2*thetas**2/pi)
for i in range(128):
As = 16.*c/pi * dot(Ks, deltas*taus)
Z = 4.*pi*sum(deltas*exp(-As))
psis = 1./Z * exp(-As)
taus = psis.copy()
S = 4.*pi*sum(deltas*psis*(3.*cos(thetas)**2-1.)/2.)
theint = dot(deltas*psis, dot(Ks, deltas*psis))
p = c + 32.*c**2*theint
pder = 1. + 64.*c*theint
f = log(c) - 1. + 4.*pi*sum(deltas*psis*log(psis)) + 32.*c*theint
mu = f + p/c
muder = pder/c
return p, pder, mu, muder
p = c + c**2
pder = 1. + 2.*c
mu = log(c/(4.*pi)) + 2.*c
muder = 1./c + 2.
return p, pder, mu, muder
def totalfunc(ct, c):
pt, pdert, mut, mudert = thefunc(ct)
p, pder, mu, muder = isofunc(c)
F = -array([pt - p, mut - mu])
J = array([[pdert, -pder], [mudert, -muder]])
dc = solve(J, F)
return ct+dc, c+dc
ct = 4.
c = 3.
for i in range(100):
ct, c = totalfunc(ct, c)
thefunc computes the self consistent solution for $\psi$ returning the pressure, the chemical potential and their derivatives.
isofunc does the same for the isotropic phase.
totalfunc does a Newton-Raphson iteration. Sorry for the names of the functions. The variables should be self explanatory, though.