ベクトル場を書いて分岐現象を見たい

ホップ分岐やサドル・ノード分岐など分岐現象を可視化する方法として方程式のベクトル場を見ることがあります．
以下はサンプルコードです．
C++コード

    
    
      #define _USE_MATH_DEFINES
      #include <math.h>
      #include <fstream>
      #include <iostream>
      #include <thread>
      #include <vector>
      const double Delta=0;
      const double deltaOmega=1;
      const double Omega1=1;
      const double Omega2=Omega1+2*deltaOmega;
      const double strength=1;
      const double p=1;
      //x軸(極座標なら動径)方向の微分方程式
      float funcx(float x, float y){
       return (-Delta + strength / 2.0*(1 - p)*(1 - x*x))*x + strength*p*(x*(1 - x*x) / 2.0)*cos(y);
       }
      //y軸(極座標なら偏角)方向の微分方程式
      float funcy(float x, float y){
       return fabs(Omega2 - Omega1) - strength*p*(1 + x*x)*sin(y);
      }
      //ベクトル場の出力
      void vector_field(){
       double beginx=0; //x軸の最小値(極座標なら動径)
       double endx=1; //x軸の最大値
       double beginy=0; //y軸の最小値(極座標なら偏角)
       double endy = 2.0*M_PI; //y軸の最大値
       double step=0.1; //ベクトル場のステップ幅
       float x, y, deltax, deltay, norm; //変数
       char filename[100];
       sprintf(filename, "vector_field.dat");
       ofstream output(filename); //ファイル出力
       for (x = (beginx-step); x <= (endx+step); x = x + step){
        for (y = beginy; y <= endy; y = y + step){
         norm = sqrt(funcx(x, y) * funcx(x, y) + funcy(x, y)*funcy(x, y));
         deltax = funcx(x, y) / (1.0 + norm)*step;
         deltay = funcy(x, y) / (1.0 + norm)*step;
         output << y << " " << x << " " << deltay/step << " " << deltax/step << endl;
        }
       }
    output.close();
    }
    int main(){
     vector_field();
     return 0
    }
    

Gnuplotコード

    
      set isosamples 200,300
      set term pngcairo enhanced size 640,480
      set polar
      set grid polar
      set size square
      set xrange [-1.1:1.1]
      set yrange [-1.1:1.1]
      set key outside top center reverse Left
      set o sprintf("vector_field.3f.png")
      plot \
      sprintf("vector_field-p%.3f-w%.3f-d%.3f.dat",p,i*a,d) with vector,
    

参考サイト

gnuplotで二次元のベクトル場を綺麗に書く。ただしベクトル場はファイルから読み込む。poissondd.wordpress.com