

function [s,t,y,dy,er,info]=shooting(F,a,b,ya,yb,s0=0.0,tolr=1e-6,tola=1e-5)
#solves Boundary value problem:
# y''=F(t,y,y') y(a)=ya y(b)=yb
#Input:
#F=handle to function:
#         function y=F(t,y,dy)
#a,b - ends of the interval (a<b)
#ya,yb - boundary values at the ends
#s0- starting value 
#tola,tolr - tolerances if |yb-y(t;s)|>tola+tolr*yb then there is a warning and info=1
#Output:
# s- the right value of s i.e .the solution y(t;s) with y'(a)=s solves BVP
# t - the mesh generated by linspace(a,b))
# y- vector with solution (1 column matrix)  on the mesh  t
# dy- vector with derivative of the solution on the mesh  t
#er - |yb-computed solution at b|
#info - 0 if |y(b;s)-yb|<=tola+tolr*|yb|
#       1 otherwise 
#Example:
#
#  F=@(t,y,dy) y*y;
#  [s1,t,y1,dy1,er1,info]=shooting(F,0,1,1,-4, -6);s1
#  [s2,t,y2,dy2,er2,info]=shooting(F,0,1,1,-4,-30);s2
#  plot(t,y1,t,y2); 
#two different solutions!
global aa bb yya yyb FF
aa=a;bb=b;yya=ya;yyb=yb;FF=F;

[s,ys]=fsolve(@G,s0);

er=abs(ys);
if(er>tolr*yb+tola)
 printf("Warning. |y(b;s)|=%8.3g > given tolerances!\n",er);
 info=1;
 return;
else 
 info=0;
 t=linspace(a,b);
 x0=[ya;s];
 y=lsode(@vecfield,x0,t);
 dy=y(:,2);
 y=y(:,1);
endif
endfunction


function y=vecfield(x,t)
 global FF

 y=zeros(size(x));
 y(1)=x(2);
 y(2)=FF(t,x(1),x(2));
endfunction

function [y]=G(s)
global aa bb yya yyb
 t=[aa;bb];
 x0=[yya;s];
 yy=lsode(@vecfield,x0,t);
 y=yy(2,1)-yyb;
endfunction

