ニュートン法による陰関数描画の並列計算

前提・実現したいこと

陰関数描画をニュートン法を用いて行い、それを並列計算させるという課題に取り組んでいます。実際にプログラムを書いて実行してみたところ、描画をするための"ans_0"～"ans_3"までは作られて図形をプロットすることはできました。しかし、実行時間が表示されず実行後に次のようなメッセージが出てしまいます。
一体このメッセージは何を表していて、どう改善すべきなのでしょうか。
よろしくお願いいたします。
なお、この問題で用いている陰関数はハートを書くもので
(x^2+y^2-1)^3-(x^2)*(y^2)=0
です。

発生している問題・エラーメッセージ

*** glibc detected *** ./pfigure: double free or corruption (out): x00007f7a3c0008c0 ***
中止

該当のソースコード

#include<stdio.h>
#include<math.h>
#include<omp.h>
#define DIV 100
#define MAXITER 10
#define ZERO 0.0
#define ONE 1.0
#define TWO 2.0
#define THREE 3.0
#define EPS 1.0E-15
double fxy(double x,double y){
  return (x*x+y*y-ONE)*(x*x+y*y-ONE)*(x*x+y*y-ONE)-x*x*y*y*y;
}
double dfx(double x,double y){
  return THREE*(x*x+y*y-ONE)*(x*x+y*y-ONE)*(TWO*x)-TWO*y*y*y*x;
}
double dfy(double x,double y){
  return THREE*(x*x+y*y-ONE)*(x*x+y*y-ONE)*(TWO*y)-THREE*x*x*y*y;
}
double newtonx(double x,double y){
  int i;
  double x2;
  for(i=0;i<MAXITER;i++){
    x2=x-fxy(x,y)/dfx(x,y);
    if(fabs(x2-x)<=EPS){
      return x2;
    }
    if(fabs(x2)>TWO){
      return x2;
    }
    x=x2;
  }
  return -THREE;
}
double newtony(double x,double y){
  int i;
  double y2;
  for(i=0;i<MAXITER;i++){
    y2=y-fxy(x,y)/dfy(x,y);
    if(fabs(y2-y)<=EPS){
      return y2;
    }
    if(fabs(y2)>TWO){
      return y2;
    }
    y=y2;
  }
  return -THREE;
}
int main(){
  int i,j;
  double unit,x,y;
  FILE *fp;
  char filename[10];
  unit=TWO/DIV;
  double t1,t2;

  t1=omp_get_wtime();

#pragma omp parallel private(x,y,i,j,filename)
  {
    sprintf(filename,"ans_%d",omp_get_thread_num());
    fp=fopen(filename,"w");

#pragma omp for
    for(i=-DIV;i<=DIV;i++){
      for(j=-DIV;j<=DIV;j++){
    x=newtonx(unit*i,unit*j);
    if(fabs(x)<=TWO)
      fprintf(fp,"%f %f\n",x,unit*j);
    y=newtony(unit*i,unit*j);
    if(fabs(y)<=TWO)
      fprintf(fp,"%f %f\n",unit*i,y);
      }
    }

    fclose(fp);
  }
  t2=omp_get_wtime();
  printf("TIME=%10.5g\n",t2-t1);

  return 0;
}