There are many different methods. Which one might be best for your particular application would really depend on the form of the equations.  Are these initial value or boundary value problems? Integro-differential equations? Integral equations? Is this anything like what's described in the paper, "Analytical Solution of a Walters’ Liquid B Flow over a Linear Stretching Sheet in a Porous Medium" by Mahabaleswar and Saha? That would make it a boundary layer problem. If so, there are several practical methods of solving boundary layer problems, starting with something as simple as the integral method.

Dear Dudley J Benton,

 The resultant differential equation is of the for $L(u)+N(u)=0$ ,where $L(u)=fxx +cfx+c1f.

and $N(u)=ffx  +f2fxxx+2ffx+fxxxx$ with two initial condition and two boundary condition.$f(0)=0$,$fx(0)=2$ , $f(b)=0$ ,$fx(b)=0$.Where b is the left boundary .and all ci are constant.Thank you for your interest.

Wow! A real nightmare. I would judiciously start analyzing the strength of influence for each term and to try like hell to remove as many as possible, then I would linearize the equations around some point of interest, finally, solve the set of equations near that point of interest. You are highly unlikely to get a general solution for such an intractable set of equations, plus, numerical techniques should be applied with extreme caution, for numeric error may corrupt any solution generated. 

With such high-order derivatives I would rather appeal to spectral methods. As pointed out by Luisiana, you can relief the burden by first eliminating terms that contribute weakly to the solution. Alternatively, you could search for some sort of simplification that allows the derivation of a "low-order" eigenfunction.

Hammad Khalil, please correct rhs to the eq. 

Umar Khan, it's naturally to try ordinary power series instead you propose, like in Uni we done, or using, say, Lie group generator, due to their high formalism. 

But, Hammad Khalil, this series should not be absolutely convergent, though with fine local asymptotics.  

Thank you , D. Kovriguine,  The RHS of the equation is L(f) and N(f,fx,fxx,fxxx,fxxxx)

Thanks, Hammad. Four order, no singularity - nice set. I've got now a similar 3D Hamiltonian problem over the solitonica, it seems of integrable case, but a series representation is hard to prove stability, since a solution with a high gradient - conventional convergency, though this one likely takes place, if follow numerical approach. Best wishes

It's not "highly nonlinear" its is very slightly nonlinear.

It is linear in all the derivatives so it is classified as a semilinear pde.

The functions multiplying the derivatives are all f or f^2 i.e. very simple.

The highest overall power of f and the derivatives is 3 so if you use a fourier method you only triple the wavelength. High powers of $f$ or functions with infinite series such as sin(f) would be much worse. Then you could have sqrt(f), i.e. bracnh points or even something non-analytic such as |f|

Given this form numerical solution using fourier series i.e. pseudo-spectral should be straightforward, or use finite differences for simplicity but less accuracy. In either case since the problem is elliptic you will need to calculate the residuals and embed the probem in an iterative solver. For efficiency you can directly solve for the linear parts.

Thank You  J .N Mcelwaine,  

Actually the problem is totally new for me and I do not understand from where I start .So if you gave me a starting point then it will be helpfull for me. I am trying adomain decomposition to solve the problem but upto know  I have not obtained satisfactory results.

The way you have formulated the problem with N and L is the Adomian method.

You should see from Wikipedia however that though this is a useful approach mathematically its not useful physically because general boundary conditions for elliptic problems will not be compatible.  

Dear Jimmy, I've reduced the order of original eq in unity. Would your pirouette be next? 

Eq:

${\frac {d^{2}}{d{x}^{2}}}f \left( x \right) +c{\frac {d}{dx}}f \left(

x \right) +c_{{1}}f \left( x \right) +2\,f \left( x \right) {\frac {d}

{dx}}f \left( x \right) + \left( f \left( x \right) \right) ^{2}{

\frac {d^{3}}{d{x}^{3}}}f \left( x \right) +{\frac {d^{4}}{d{x}^{4}}}f

\left( x \right)$

Reduce:

$\left( {\frac {d^{3}}{d{a}^{3}}}b \left( a \right) \right) \left( b

\left( a \right) \right) ^{3}+ \left( {\frac {d}{da}}b \left( a

\right) \right) b \left( a \right) +cb \left( a \right) +c_{{1}}a+2

\,ab \left( a \right) +{a}^{2} \left( {\frac {d^{2}}{d{a}^{2}}}b

\left( a \right) \right) \left( b \left( a \right) \right) ^{2}+{a

}^{2} \left( {\frac {d}{da}}b \left( a \right) \right) ^{2}b \left( a

\right) +4\, \left( {\frac {d^{2}}{d{a}^{2}}}b \left( a \right)

\right) \left( b \left( a \right) \right) ^{2}{\frac {d}{da}}b

\left( a \right) + \left( {\frac {d}{da}}b \left( a \right) \right)

^{3}b \left( a \right) =0$,

where

$\left\{ a=f \left( x \right) ,b \left( a \right) ={\frac {d}{dx}}f

\left( x \right) \right\} $.

Dear Dmitry Kovriguine,

                                      The idea of reduction of the order of derivative is amazing , but I  do not understand that the reduction of order will ease the way to approximate solution, because as you see that before reduction only three terms are nonlinear but after reduction we have I think six terms which are nonlinear.

In general reducing the order is good because the resulting matrices are better conditioned so can solve them more easily and with greater accuracy. It depends really on whether you have an IVP or BVP. If you have an IVP then you can just use a standard ODE solver otherwise you probably need to stick with the 4th order equation and solve iteratively

Dear colleagues, we always should suspect to reduce the problem to a simplest case. For instance, who can say that third-order algebraic curves are not fiction, if not trying to reduce them, using some transform, to those of the second order; ellipses, parabolas etc. But if these, the thirds, exist then we can study them. Who says that we cannot reduce the order in original equation one more time?

If I may... higher order derivatives in time or space are interesting only on small time or spatial intervals, otherwise, their influence forces a far too complex global solution space; hence, it is usually impossible to generate a global sol'n for such equations. The question of reducing the degrees of freedom is more a question of "what exactly is the question?", eh? In other words, what is of interest to be studied... It may be that the lower order derivatives need to be thrown out or vice versa. 

The main quiestion is it possible to approximate the solution $f(x)$  at least on some finite domain.

Hello, Hammady and Jimmy. I've got three first approximations (in the memory of the Great October Red Rev.):

$\left\{ [f \left( x \right) =\sin \left( 4.730040745\,x \right) -

1.017809411\,\cos \left( 4.730040745\,x \right) + 0.008904705330\,{

{\rm e}^{ 4.730040745\,x}}+ 1.008904706\,{{\rm e}^{- 4.730040745\,x}},

\left\{ a=- 4516137849000.0- 9249428545.0\,b \right\} ],[f \left( x

\right) =\sin \left( 10.99560784\,x \right) - 1.000033551\,\cos

\left( 10.99560784\,x \right) + 0.00001677550011\,{{\rm e}^{

10.99560784\,x}}+ 1.000016776\,{{\rm e}^{- 10.99560784\,x}}, \left\{

a= 8916543595000.0+ 614140937.7\,b \right\} ],[f \left( x \right) =

\sin \left( 17.27875966\,x \right) - 1.000000063\,\cos \left(

17.27875966\,x \right) + 0.00000003132781030\,{{\rm e}^{ 17.27875966

\,x}}+ 1.000000031\,{{\rm e}^{- 17.27875966\,x}}, \left\{ a=-

175185138700000.0- 1971220432.0\,b \right\} ] \right\}$

Hello,Khalil. I have used the GA with orthogonal polinomial for this kind of problem instead of linearization's Newton Method, since this one use to diverge.

I have the code in Matlab ready to work, if you wish.

All the best,

Osvaldo