#include <stdio.h>
#include <math.h>
#define  EPS 1.0e-6
double f(double x, double y){return x*2/3 + (x*x - y*y)/(3*((x*x+y*y)*(x*x+y*y)));}
double g(double x, double y){return y*2/3 + 2*x*y/(3*((x*x+y*y)*(x*x+y*y)));}

void newton_hukuso(double x, double y, double eps, double *result_x, double *result_y)
{
    int i=0;
    double nx, ny;

    while(i<8) {
        i++;
        nx = f(x, y);
        ny = g(x, y);
        // 収束条件を満たせばループ終了
        if(fabs(nx - x)<eps) break;
        x = nx;
        if(fabs(ny - y)<eps) break;
        y = ny;
    }
    *result_x = x;
    *result_y = y;
}
int main()
{
    double x = -2, y = -2, result_x, result_y;
    FILE *fp;
    fp = fopen("plot2.dat", "w");
    for(x=-2; x<=2; x+= 4/5){
      for(y=-2; y<=2; y+= 4/5){
        if(x==0 && y==0){
          break;
        }
        // ニュートン法(初期値, 収束条件)
        newton_hukuso(x, y, EPS, &result_x, &result_y);
        // 解を書き込む
        if(sqrt((result_x-1)*(result_x-1) + result_y*result_y) <= EPS){
          fprintf(fp, "%.13lf %.13lf\n", result_x, result_y);
        }
      }
    }
    fclose(fp);
    return 0;
}