@Lewis: I noticed that the integrands in both of your integrals become infinite at the upper limit of integration, and this can result in some inaccurate numerical results from the integration process. However, with the appropriate change of variables these integrals can be converted to complete elliptic integrals of the first and second kinds for which Matlab has functions that you can call on directly and thereby avoid using numerical integration at each step of the iteration.
Here is an outline of that change of variable. You are working in a realm where lambda and c must both be negative with lambda having the greater absolute value. Define A = -c and B = -lambda, so that A and B are positive quantities with A < B. Make the following change of variable:
y = sqrt(B-A)*sin(t)
dy = sqrt(B-A)*cos(t)*dt
(After all the smoke clears away) the first of your integrals becomes the integral with respect to t from 0 to pi/2 of
which can be evaluated in terms of a complete elliptic integral of the first kind. It is easy to show that your second integral will reduce to a combination of complete elliptic integrals of the first and second kinds.
Hence in computing your objective function for 'fsolve' you need not call on Matlab's 'integral' function at all but only on the complete elliptic integral functions.