% Usage: 
% y = linshoot(p,q,f,a,b,ya,yb,M,n) or
% [y t] = linshoot(p,q,f,a,b,ya,yb,M,n) or
% [y t s] = linshoot(p,q,f,a,b,ya,yb,M,n)
% Shooting method for linear boundary value problems of form y"+py'+qy=f
% using a fixed step-size method given by M
%
% Input:
% p,q,f - Matlab inline functions p(t), q(t), f(t)
% a,b - interval
% ya, yb - boundary conditions
% M - fixed step-size solver
% n - number of steps
%
% Output:
% y - computed solution
% t - time steps
% s - value of y'(a)
%
% Examples:
% [y t]=linshoot(inline('0'),inline('1'),inline('0'),0,1,1,1,@rk4,10);

function [y t s] = linshoot(p,q,f,a,b,ya,yb,m,n)
% compute v, first expressing it as a solution to first order eq.
fu = @(t,u) [u(2), f(t) - p(t) * u(2) - q(t) * u(1)];
ua = [ya 0]; % initial condition for u: u(a)=ya, u'(a)=0
[u t] = m(fu,a,b,ua,n);
% compute u
fv = @(t,v) [v(2), - p(t) * v(2) - q(t) * v(1)];
va = [0 1]; % initial condition for v: v(a)=0, v'(a)=1
[v t] = m(fv,a,b,va,n);
% check if v(b) is too small
if abs(v(n+1,1)) < 1e-13
    if u(n+1,1) == yb
        display('The problem may have infinitely many solutions.');
        display('Warning: The provided solution might be inaccurate!');
        y = u(:,1);
        s = 0;
    else
        display('The problem may have no solution.');
        display('Warning: The provided solution is useless.');
        y = zeros(length(t));
        s = 0;
    end;
else
    s = (yb - u(n+1,1)) / v(n+1,1);
    y = u(:,1) + s * v(:,1);
end;