Fourth order differential equation using bvp4c gives singular jacobian error

I am trying to solve a fourth order differential equation of the following form, using bvp4c: d4y/dz4=-p(z,y) subject to the boundary conditions: y'''(LPile)=0; y'''(0)=F/ModE*Inertia; y''(LPile)=0; y''(0)=M/ModE*Inertia;

This is the code I've come up to.

dzdy=@(z,y)[y(2);y(3);y(4);(-1/(ModE*Inertia))*gamma*z*D*2.4*((qc/(gamma*z))^0.67)*((z/D)^0.75)*(1-exp((-6.2)*((z/D)^(-1.2))*((y(1)/D)^0.89)))]; res=@(ya,yb)[ya(4)-(F/ModE/Inertia);yb(4);ya(3)-(M/ModE/Inertia);yb(3)];

solinit = bvpinit(linspace(0,LPile,2*LPile),[0.1;0.1;0.1;0.1]);

Fsol = bvp4c(dzdy, res, solinit);

 But it gives me the following error:

??? Error using ==> bvp4c at 252 Unable to solve the collocation equations -- a singular Jacobian encountered

Error in ==> LatestCode20150721 at 28 Fsol = bvp4c(dzdy, res, solinit);

z appears 4 times in the equation of p(z,y), which I guess is the problem, as the code runs fine when I use less complicated functions for p(z,y). Any ideas how to fix this?

In reality qc is also a function of z; which I just kept out for simplicity. The input data I am using are as follows: LPile=35; ModE=2*(10^8); gamma=20; D=4; Inertia=10.32; F=-10000; M=40000; qc=5;

[Torsten]

Maybe y(1) becomes negative in the course of the iteration such that (y(1)/D)^0.89 produces a complex result. Just a guess ...

Best wishes

Torsten.