

function [s,t,y,dy,er,info]=linshooting(d,c,f,a,b,ya=0,yb=0,tolr=1e-4,tola=1e-6)
#we are solving BVP
#-y''+d(x)y'+c(x)y=f with y(a)=ya;y(b)=yb
#Input:
#f,d,c=handle to functions:
#         function y=f(x)
#         function y=d(x)
#         function y=c(x)
#a,b - ends of the interval (a<b)
#ya,yb - boundary values at the ends
#
#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 and its derivatives (2 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:
#d=@(x) 0;f=@(x) x; c=@(x) 1; a=0;b=1; #y''=y+f y(0)=y(1)=0
#[s,t,y]=linshooting(d,c,f=@(x) x,a,b);plot(t,y,";y''=y+x,y(0)=y(1)=0;");
#
#d=f=@(x) 0; c=@(x) 1; a=0;b=1; #y''=y y(0)=y(1)=1
#[s,t,y,dy,er,info]=linshooting(d,c,f,a,b,ya=1,yb=1);plot(t,y,";y''=y,y(0)=y(1)=1;");
#test with  b = 5,10,15,20, is the error growing? Why?

global cc dd ff
dd=d;cc=c;ff=f;

tol=tola+tolr*yb;
tt=[a;b]; 
x0=[ya;0]; #s_1=0
y1=lsode(@vf,x0,tt);
y1=y1(2,1);
x0=[ya;1]; #s_2=1
y2=lsode(@vf,x0,tt);
y2=y2(2,1);
 
if (abs(y2-y1)> 1e-14*norm([y1,y2],"inf"))
   s=(yb-y1)/(y2-y1);
else
   printf("\nWarning: |y(%5.2g;1)-y(%5.2g;0)| is almost zero. terminating\n\n",b,b);
   return;
endif
if(nargout >1)
 t=linspace(a,b);
 x0=[ya;s]; #s
 y=lsode(@vf,x0,t);
 er=abs(y(100,1)-yb);
 dy=y(:,2);
 y=y(:,1);

 if(er<=tola+tolr*yb)
   info=0;
 else 
   info=1;
   printf("\nWarning. Computed s=%g is wrong, as |y(b)-yb|=%g which is above given   tolerances\n\n",s,er);
 endif
endif #nargout>1
endfunction
 

function y=vf(x,t)
 global cc dd ff
 #x_1=y;x_2=y' vector field
 #f(x,t)=[x(2);d(t)*x(2)+c(t)*x(1)-f(t)]
  y=zeros(size(x));  
  y(1)=x(2);
  y(2)=dd(t)*x(2)+cc(t)*x(1)-ff(t);
end
