thanks

it solved with following code :

dydx=@(x,y)[y(2) ; y(1)^2+y(5)*y(2)-y(3)^2 ; y(4) ; 2*y(1)*y(3)+y(5)*y(4) ; -2*y(1) ];

BC=@(ya,yb)[ya(1) ; yb(1) ; ya(3)-1 ; yb(3) ; ya(5)];

yint=@(x)[0.4 ; 0 ; 0.2 ; 0 ; .6 ];

solint=bvpinit(linspace(0,20,2000),yint);

S=bvp4c(dydx,BC,solint);

for i=1:4

S=bvp4c(dydx,BC,S);

end

In cases where an analytical solution is not obtainable, a very simple technique is the finite-difference method. See Lecture 9 here:

http://emlab.utep.edu/ee5390cem.htm.

You missed a negative sign defining y(4). It should be 2*y(1)*y(3)-y(5)*y(4).