Fist of all the integration term has to be reduced by taking the differentiation in the entire equation w.r.t. z once. Finally, we obtain a third order differential equation which is free from  the integration term.

This equation can be solved using iteration. Firstly, by solving without integral term we can obtain  solution for  first iteration.  For second iteration solution we can use the first iteration solution for evaluation the integral term in this equation and so on.  

Problems like this (integro-differential equations) often arise in EE/circuit analysis. The conventional treatment is to apply a suitable integral transform (e.g. Laplace transform), whereby the integral disappears, then solve the resulting equation in "transformed" space, and finally use the inverse transform to recover the physical solution. I suspect this will be tricky since your equation is non-linear. I think a good place to start would be the text by Lakshmikantham and Rao (1995) Theory of Integro-Differential Equations or something similar.

Use Method of Variation of Parameter to solve second order differential equation.

OR Use by Complementary Solution and Particular Integral Method.

One complication in using the  finite difference method for solving the boundary value problem is that the integral (say approximated by Simpson's rule) makes the matrix to be inverted will become non-sparse (since the "stencil" at every point spans the whole domain now)  rather than three points. The other  issue is that the nonlinear term BP^3....going by what you said I am assuming you have found a way to handle the BP^3 term iteratively. The integral being of P and [not any f(P)] does not seem to add nonlinearity to the problem. So, hopefully u can find an iterative matrix inversion scheme that accomodates both issues.

If you want to get the numerical solution you have not any problem.

I think you can follow the following steps:

The "variation of parameters" method was suggested above, but this is applicable only for linear ODEs (this one is non-linear).

Thank you for all professors.  Dear Dr Zakaria, Thank you for your contribution. I have looked at you answer. The parameter P is a a function of Z. You have changed the 2nd ODE to 1st ODE  and this is nice. But the integration still in the equation which is P(z)dz. even in the numerical method how to fine the integration P(z)dz if I am seeking P(z).

Dear Dr. Ahmad 

I ask Allah to save you and your family.

As I said before the integral term is a constant value C0, but we seek to evaluate it  by using an appropriate  numerical method, so I suggest the following procedure:

Note :

The derivation of the exact solution not just reduction the order of ODE.

The integration of two order equations into integral equations by using the fixed point theorem, then can solve

Dear Dr. Ahmad,

Solving this ODE seems easy. If you have Matlab Or maple I can provide you a subroutine to solve the equation. What about boundary conditions?

Yes Dear Dr. Mohammad. 

Thank you very much for your answer. the boundary condition is (dp/dz)=+-p/delta at -+L/2 where delta is constant.