You should specify more details about your problem.

The general principle of preparing bifurcation diagrams (bifurcation parameter vs. state variable) is as follows:

On x-axis put a bifurcation parameter and on y-axis the values of the state variable of interest, you have obtain while doing N Poincaré sections (e.g. N=100) .

There are various matlab codes to generate bifurcation diagrams at Mathworks File Exchange.

For example, this one generates bifurcation diagrams for 1-D maps:

https://www.mathworks.com/matlabcentral/fileexchange/26839-1d-bifurcation-plot?requestedDomain=www.mathworks.com

I believe you must use a more-stable integrator, so you obtain non-spurious solutions. Try for example, implicit o semi-implicit Runge Kutta. See:

https://www.academia.edu/11949052/Implicit_Runge-Kutta_Algorithm_Using_Newton-Raphson_Method

https://www.academia.edu/11948884/Lobatto_Implicit_Sixth_Order_Runge-Kutta_Method_For_Solving_Ordinary_Differential_Equations_with_Stepsize_Control

Code in fortran RKI36 here in my RG publications.

Dear Michał Marszal,

Thanks for your response. I plotted the bifurcation diagram and enclose it here. The behavior is right but it is too thick and some extra information. This bifurcation diagram showing here is for investigations of the Driven, Damped, Nonlinear Pendulum.By increasing the N, figure become thicker but loose some information. Do you have any idea?

Bests.

Dear Tamas Kalmar-Nagy 

Thanks for your response.

Do you please have any idea to use this code in Lorenz system or Driven Pendulum? The Lorenz system have 3 differential equations that attached below. 

Bests.

Dear Andrés Granados

Thanks for your response.

The articles sound great and I am going to read them.

bests.

Hi Dear Abd AlRahman 

Thanks for your nice response.I think you are very close to answer. Would you mind if I have your Matlab code? I wrote it in rk4 and it works right. but i have problem in this bifurcation diagram.

thanks again.

Bests

Your bifurcation diagram looks strange. Are you marking each state variable with different colors? You usually mark single state variable in y-axis (e.g. chosen variable or norm of the state vector) How do you make Poincaré sections (stroboscopic map or while one of the state variables is crossing zero)?  Is your system of ODEs autonomous or non-autonomous?