I suppose it's the biharmonic equation for some beam problem? Shooting is probably not the best choice. A finite difference or finite element method is in all likelihood preferable.

If it is a beam problem, and you can "decompose" the problem into a moment equation M" = q, and a deflection equation u" = M/(EI), this formulation is much preferable. Discretizing the fourth order operator directly is more sensitive to perturbations and roundoff, than if you use the "decomposition" above, where you only have to discretize second order operators (but instead solving two of them).

If it's this problem, you have lots of finite difference (FDM) and finite element methods (FEM) to choose from. For example, you could use a 4th order convergent FDM on a uniform grid and get very high accuracy at little cost. If you need a nonuniform grid, it is better to use FEM, as this is more flexible with nonuniform grids.  

Dear prof. Gustaf,

   Thank you for your insightful reply, I attached a summarized description of the problem which arises in peristaltic flow dynamics.

I wait to listen your helpful suggestions. 

With best regards

Yours, Khaled G. Mohamed

Ok, a very different application, but nevertheless it could be factorized into two second order problems as I suspected. Then, by all means, avoid solving the forth order problem directly, and instead focus on solving the two second order problems. Because the first one is integrable by analytic techniques (solution provided in your problem description) there is only one second order problem to solve, really, so everything is pretty standard.

The difficulty here lies in the nonlinearity, because no matter how you discretize the problem, you will have to use a nonlinear  (Newton) solver to iterate. That takes a good starting guess of the solution. Also, it is not clear whether you can do this problem on a uniform grid, but that's where I would start. (Today it is certainly affordable to solve problems also on a very fine grid, say having up 10,000 mesh points.)

The description's concern over how one can solve this problem for unknown integration constants lambda1 and lambda2 seems to be the least concern of all.

However, don't use shooting -- it is likely to get you into more trouble than the problem warrants. Start out using a straightforward 2nd order accurate finite difference two-point boundary value problem solver. I don't know about your preferences, but working in Matlab is usually quick. Your greatest trouble will be to get the boundary conditions right to second order accuracy, and get the Newton iteration (due to the nonlinearity) to converge.

Best of luck!

Oops, made a mistake here. The problem is not nonlinear as y is the independent variable. So the problem is very straightforward. Just use a standard finite difference method, only pay attention to getting the Neumann boundary conditions right.

Otherwise the above remarks are valid. Don't use shooting (because you will encounter stability problems and will have to iterate), and don't bother with the fourth order problem. You have an awful lot of analytical information available, do don't let it go to waste.  

Gustaf has provided some valuable suggestions and good pointers. An yes, going the finite difference method is good suggestion especially if you are  comfortable formulating the finite difference code. As you point out the  problem (peristaltic flow) you want to solve is likely to be nonlinear so that means you are going to have to incorporate a Newton method into your solution technique. But you are dealing with a set of nonlinear ODEs so  getting good initial guesses for the Newton method is then the key.  Grid refinement may be  another important issue. Things become more interesting and challenging if you want to see how the solution evolves as you change a parameter- so-called continuation analysis. This can be done by incorporating a continuation method with your BVP,  but in some cases solving the PBP using a shooting method may be more efficient.  For example, you do not have to worry about automatic grid refinement which is built into the MatLab of Mathematica ODE solvers.  If you would like to see how this is done with a coupled set of two nonlinear BVP, check out  my report "Computing Hysteresis Effects in Plane Poiseuille Flow with Viscous Heating". All these issues are addressed in this report. Good luck.

There are various ways of solving your equation as has already been touched upon. However, I am a great advocate for the shooting method, and unless your equation is stiff, then it will be very readily amenable to solution by the shooting method (assuming that one has a sufficiently powerful code, i.e. one which will solve any order of system and any number of unknown initial conditions).

Your second order version with two unknown constants may be solved this way too. If you create four variables in this way: v_1=psi, v_2=psi', v_3=lambda_1, v_4=lambda_2, then four first order equations may be written down. The third and fourth of these are (v_3)'=0 and (v_4)'=0, since both these 'variables' are constant. The two unknown initial conditions are v_3(0) and v_4(0), while the two conditions to be satisfied at the far end of the interval are the values of v_1 and v_2. This should work very nicely, and this 'trick' is something I have used for decades.

A very recent paper of my own where such a procedure is the following: D.A.S.Rees, A.P.Bassom, G.Genç (2014) Weakly nonlinear convection in a porous layer with multiple horizontal partitions. Transport in Porous Media 103 437-448. This is a convective instability analysis for a system of thermally coupled porous layers. We used multiple shooting, so that the final values at the end of one sublayer are related to the initial values at the start of the next sublayer, i.e. all the sublayers were computed in parallel rather than in sequence. A critical Rayleigh number is computed in the same way as one of your lambda values could be. The neutral curve was minimised, hence the wavenumber was found in the same way as the second of your lambda values. Weakly nonlinear theory meant that many further systems of ODEs needed to be solved. So for a system comprising of N layers, there are 36N+5 equations to be solved (the 5 corresponding to constants like your lambdas) with 36N-9 unknown initial conditions. We managed to compute for N=100 layers at the most (i.e. 3605 equations with 3591 unknown initial conditions), but this was because of the amount of time taken to solved such a system.

My conclusion is for you not to give up on the shooting method. It is incredibly versatile, particularly if you have written your own black-box code, and can be used to solve some very very large systems of equations which a direct solution using standard finite differences would take a long time to set up. You do need a general code, though, one which isn't devoted to systems of a fixed order. As long as your system isn't stiff, then it ought to work well. My feeling is that there should be no essential difference in difficulty in solving your 4th order version or your 2nd order version with two constants.

Dear Profs. Gustaf, Brian and D. Andrew,

    Thank you for your helpful reply, but what is your opinion about using the following numerical method in which the 4th order boundary value problem may be directly solved.

with my best regards

Viswanadham et al. Numerical Solutions of Fourth Order Boundary Value Problems by Galerkin Method with Quintic B-splines, International Journal of Nonlinear Science

Vol.10(2010) No.2,pp.222-230  

Included as attachment.

Much depends on the library codes or the bespoke codes which you have available to you for your use. I have my own shooting method code which I created 28 years ago and which I have used for a huge variety of different ODE problems. Other direct methods are also excellent, especially for linear equations such as yours.  A 4th order system is pretty small and therefore there is little to choose between the various methods. That said, if you were to use a shooting method, then the basic solver could be a 4th order Runge Kutta scheme, or better, whereas most direct solvers using finite differences are often only of second order accuracy only. I don't know how accurate quintic B-splines are, but they may well be very good, and equal to a 4th order Runge-Kutta scheme. You will need to answer the question: which method will be the fastest to implement? Even if the method you choose to use is of second order accuracy, you will be able to use Richardson's extrapolation to improve the order of accuracy. Successive interval-halving will always indicate how many significant figures of accuracy you have.

I apologise for being a little vague about the different methods, but I have also used direct methods such as Prof Gustav has mentioned, and they too work very well indeed.  There is little to choose between these methods in my view, and therefore the speed of implementation of the method you personally choose to use (i.e. the time taken to create a correct working code) may well be the most important factor to consider.

Dear Khaled Mohamed,

I do not think, that it is a good idea to solve numerically such simple equation. But in the case of more complex equation and boundary value problem You can use tridiagonal matrix algorithm or invariant embedding methods. For 2n order with n>1 You can use matrix variant. It simple idea, R=U'/U so You have first order equation (initial value problem) for R and then another one first order equation (again initial value problem). I do not know good book for matrix variant, but it is very simple (be careful with matrix). One example of use is here

https://www.researchgate.net/publication/258521747_A_modification_of_the_invariant_imbedding_method_for_a_singular_boundary_value_problem?ev=prf_pub

This methods work very good for boundary value problem, but only in linear case (the nonlinear variant is very complex).

Good Luck

Article A modification of the invariant imbedding method for a singu...