# Solving condensate density problem in MATLAB

I want to solve for $n_{0}$ for a fixed value of $n$, lets say $n=1$ $$n= n_{0}+ \dfrac{1}{2}\int_{-1/2}^{1/2}dq\left(\dfrac{e_{q}+Un_{0}}{\hbar\omega}-1\right)$$ where $e_{q}=2[1-cos(2\pi q)]$ and $\hbar\omega=(e_{q}^{2}+2Un_{0}e_{q})^{1/2}$

Then I want to plot between $U$ and $n_{0}$. The following is the code in MATLAB.

 eq = @(q) 2*(1-cos(2*pi*q));
hq = @(q,U,n0) ((eq(q))^2+2*U*n0*(eq(q))).^(1/2);
y = @(q,U,n0) (((eq(q))+(U*n0))/hq(q,U,n0))-1;
a = -0.5;
b = 0.5;
cv =@(U,n0) n0+(0.5*v)-1;
n0 = 0.1;
options = optimset('Display', 'iter');
n = @(U) fsolve(vv,n0,0.1);
U = 0:0.1:20;
plot(U,n(U))


edit 2 code corrected edit 1.

Edit 3

But I still facing the following errors:

Error using ^ Inputs must be a scalar and a square matrix. To compute elementwise POWER, use POWER (.^) instead.

Error in @(q,n0)((eq(q))^2+2*Un0(eq(q))).^0.5

Error in @(q,n0)(((eq(q))+(U*n0))/hq(q,n0))-1

Error in @(q)y(q,n0)

Error in quad (line 72) y = f(x, varargin{:});

• Welcome to SciComp! As currently posed, your question is unclear. As WolfgangBangerth points out, it's not clear what you've tried. It looks like you're relying on quad to calculate the integral, and then passing that integral as a function to fsolve. The algorithms used by fsolve require derivatives of the nonlinear equation being solved. Since quad uses adaptive quadrature algorithms, it would not surprise me if the finite difference derivative approximations used by fsolve were ill-behaved. You might have better luck using a non-adaptive quadrature instead. Feb 21 '15 at 5:19
• The first power in your second line needs to be eq(q).^2, not eq(q)^2. (Before you edit again: the following lines have similar problems; anytime a vector such as eq(q) appears, you need .*, ./ etc.) Feb 25 '15 at 15:32