I'm looking for a maximum accuracy quadrature formula:

$$\int_{-1} ^{1} \sqrt{\frac {1-x}{1+x}} f(x)dx = A_1f(x_1)+A_2f(x_2)+R(f)$$

I don't know exactly if it's Trapezoidal rule which has the degree of accuracy one, or Simpson's rule with three or Gauss (these 3 I have studied deeply, but might be others), also I would like to prove it using numerical analysis methods so any help in this direction will be highly appreciated.

-
I suggest Gauss-Chebyshev. See dlmf.nist.gov/3.5.E25 – Johannes Trost May 10 '14 at 22:13
Could you define "maximum accuracy" precisely? The usual requirement is that it holds for all polynomials $f$ up to certain degree. – Alexandre Eremenko May 10 '14 at 22:27
Degree of accuracy = precision of the quadrature formula. That is, as far as I have studied it, the largest positive n such that the formula is exact for x^k, k=1..n – Soreena May 10 '14 at 22:46

I assume you mean maximum degree of accuracy for polynomials, that is, you require that the formula is exact for polynomial functions of degree up to $d$, and you look for the maximum possible $d$. This is the most usual requirement, although not necessarily the best one (see for instance http://eprints.maths.ox.ac.uk/1116/1/NA-06-07.pdf).

The maximum $d$ achievable clearly depends on the number of terms $n$ in the formula $\sum_{i=1}^n w_if(x_i)$ that you are allowing for. So you are comparing apples and oranges in your examples: Simpson has a higher order than the Trapezoidal Rule, but this is compensated by the fact that it requires more function evaluations; it is just a trade-off in many cases. Note, also, that the orders that you are stating for these two rules are for quadrature without a weight function, while you have $w(x)=(1-x)^{1/2}(1+x)^{-1/2}$ added as a weight function. This changes the problem, because now "exact for polynomials" means a different thing.

It is a classical theorem, often taught in advanced numerical analysis courses, that Gaussian quadrature (for a specific weight function) has the highest possible order of accuracy among all those with a given number of terms --- see e.g. Stoer, Bulirsch, Introduction to numerical analysis Sec. 3.6.

Your weight function is a special case of the Gauss-Jacobi weight $w(x)=(1-x)^\alpha(1+x)^\beta$, so you should look for formulas giving weights and nodes for Gauss-Jacobi quadrature.

-

As Federico Poloni said, it is a special case of Gauss-Jacobi quadrature. I am in a good mood, so here is your homework. The formula is $$\int_{-1}^{1}\sqrt{\frac{1-x}{1+x}}f(x)dx= \frac{\pi(5+\sqrt{5})}{10}f\bigg(\frac{-1-\sqrt{5}}{4}\bigg)+\frac{\pi(5-\sqrt{5})}{10}f\bigg(\frac{-1+\sqrt{5}}{4}\bigg)+R(f).$$

-

Federico's answer gives you the theory. Since you only ask for two abscissae, it is easy to solve this case by first principles. Just put in the abscissae as variables and do the integral. $f(-3^{-1/2})+f(3^{-1/2})$ gives the right answer for cubic polynomials and not for quartic polynomials. I hope this isn't homework.

-