

function [s0,Y,x,erb]=linshoot(d,c,f,a,b,ya,yb,N=100)
%we solve BVP -y''+d(x)y'+c(x)y=f(x) in (a,b)
%y(a)=ya y(b)=yb by a shooting methods
%i.e. we want to find s0 s.t.
%the solution of IVP: y(t;s) 
%  -y''+d(x)y'+c(x)y=f(x) in (a,b)
%y''= d(x)y'+c(x)y-f(x) or y(1)'=y(2);y(2)'=d(x)y(2)+c(x)y(1)-f(x) 
%y(a)=ya y'(a)=s0 satisfies y(b;s0)=yb
%to compute s0 we solve numerically twice IVP for s=0,1 getting
%y(b;0) and y(b;1) then s0=(yb-y(b;0))/(y(b;1)-y(b;0))
%INPUT
%d,c,f - function handles to respective coeff. functions d(c),c(x),f(x)
%a,b - interval endswitch
%ya,yb- boundary values
%N - no of discrete poitns minus one (default 100)
%OUTPUT
%s0 value of y'(a) s.t. the solution of IVP solve BVP
%x- vector of discrete arguments (equidistant on [a,b]) 
%Y - a vector of approxmations of the BVP solution and its derivative
% at discrete equaidistant points over [a,b]
%EXAMPLE 1
%%-y''+y=0 a=0 b=1 ya=yb=1 (d=c=f)
%d=@(x) 0; f=d;c=@(x) ones(size(x));
%a=0;b=1;ya=yb=1; [s0,Y,x,erb]=linshoot(d,c,f,a,b,ya,yb,N=100);[s0,erb]
%%Example 2 -y''=y y(0)=1 y(T)=1 T=1,2,pi,pi*(1+eps) eps small
%d=@(x) 0;f=d;c=@(x)-ones(size(x));a=0;b=1;ya=yb=1; [s0,Y,x,erb]=linshoot(d,c,f,a,b,ya,yb,N=100);[s0,erb]

global dd cc ff
dd=d;cc=c;ff=f;
X=lsode(@g,[ya;0],[a,b]);
yb0=X(2,1)
X=lsode(@g,[ya;1],[a,b]);
yb1=X(2,1)
yb1-yb0
if(abs(yb1-yb0)<1e-6)
 disp("There may be no unique solution!");
end 
s0=(yb-yb0)/(yb1-yb0)
x=linspace(a,b,N+1)';
Y=lsode(@g,[ya;s0],x);
erb=Y(N+1,1)-yb;
endfunction


function z=g(y,x)
  global dd cc ff
  z=zeros(2,1);
  z(1)=y(2);
  z(2)=dd(x)*y(2)+cc(x)*y(1) -ff(x);
endfunction
