Bifurcation diagram by C++ of the non-local elastica

Question

I have tried to do a simulation using C++ of the equation of non-local elastica defined by the following system

Using the following discretization

with the boundary condition and by sitting lc=1/(6*sqrt(3)), but I did not arrive to the bifurcation diagram like in the following image,

This is the code that I have used.

    #include <iostream>
#include <math.h>
#include <fstream>
#include <cstdio>
#define PI 3.14
using namespace std;
main (){
 int i,n,m;
double a=0,b=PI,h,l,BETA=PI*PI;
 double X[500],V[500], T[500];
 cout<<"donne n=";
 cin>>n;
 l=1./(2*n*sqrt(3));
ofstream SOL("eringen.txt");
V[0]=0;
V[100]=0;
X[0]=5;
//SOL << X[0]<<"\t"<< V[0] << endl;

 for (m=1; m<=100; m++)
{
for (i=1; i<=10; i++)
 {
    X[i]=X[i-1]+(1./m)*V[i-1];
        V[i]=V[i-1]-(1./m)*BETA*(1+l*l*V[i-1]*V[i-1])*sin(X[i])*(1./(1-BETA*l*l*cos(X[i])));
        SOL << X[i-1]<<"\t"<< V[i-1] << endl;   
}
}
 FILE *fp = popen("gnuplot", "w");
fprintf(fp,"plot 'eringen.txt'  ;\  pause mouse \n");
cin.get();
 pclose(fp);
}