
function [X,t,s,er]=linshoot15(p=@(t) 0,q=@(t) -1.0,f=p,a=0,b=1,alpha=1,beta=1,N=100)
#shooting for y''+py'+qy=f
#default p=f=0 q=-1 a=0 b=1 alpha=beta=1
#with boundary values
#y(a)=alpha y(b) =beta
#input
#f,p,q - functions handle to function in the ode : y''+py'+qy=f
#a,b,- the times for bnd values
#alpha,beta - bnd values i.e. y(a)=alpha y(b)=beta
#N - no of times for which solution should be computed
#Output:
#X- the solution on [a,b] (100 points)
#t- equidistant meshpoints of [a,b] (N points)
#s - the right value of the shooting method i.e. the solution
#    with y'(a)=s solves y(b)=beta
global ff pp qq
ff=f;pp=p;qq=q;
vvf=@vf;
#a=0;b=alpha=beta=1;N=100;
#f=@(x) 0; p=f; q=@(x) -1.0;

X=lsode(vvf,[alpha;0],[a,b]);
x0=X(2,1);

X=lsode(vvf,[alpha;1],[a,b]);
x1=X(2,1)-x0;

#X(b,s)=x0+s*x1; so X(b,1)-X(b;0)=x1 X(b;0)=x0

if(abs(x1)>eps*beta);
 s=(beta-x0)/x1;
else
 printf("Warning: the solution at b=%g for y'(a)=1 is may be zero i.e. =%g\n",b,x1);
endif 

#test if s is OK
X=lsode(@vf,[alpha;s],[a,b]);
er=abs(beta-X(2,1));

t=linspace(a,b,N);
X=lsode(@vf,[alpha;s],t);
endfunction 



function y=vf(x,t)
global pp qq ff
y=zeros(2,1);

y(1)=x(2);

#y(2)=ff(t)-pp(t)*x(2)-qq(t)*x(1);
 y(2)=x(1);
end
