function w = linshoot ( coeff, a, b, n, alpha, beta )

%LINSHOOT     approximate the solution of the linear boundary value
%             problem
%
%                   u'' = p(x) u' + q(x) u + r(x)
%                   alpha(1) u(a) + alpha(2) u'(a) = alpha(3)
%                   beta(1)  u(b) + beta(2)  u'(b) = beta(3)
%
%             using the shooting method, with the solution of all 
%             initial value problems approximated using the classical
%             4th-order Runge-Kutta method
%             
%
%     calling sequences:
%             w = linshoot ( coeff, a, b, n, alpha, beta )
%             linshoot ( coeff, a, b, n, alpha, beta )
%
%     inputs:
%             coeff   string containing name of m-file defining the 
%                     functions p(x), q(x) and r(x) on the right-hand
%                     side of the differential equation; function should
%                     take a single input and return the values of p, q
%                     and r, in that order; i.e., the m-file header should
%                     be of the form
%                            [p, q, r] = coeff ( x )
%             a       left endpoint of problem domain
%             b       right endpoint of problem domain
%             n       number of uniformly-sized steps to take in marching
%                     from x = a to x = c
%             alpha   three-component vector of the coefficients which
%                     define the boundary condition at x = a
%             beta    three-component vector of the coefficients which
%                     define the boundary condition at x = b 
%
%     output:
%             w       vector of length n+1 containing the approximate
%                     values of the solution of the boundary value
%                     problem at the locations x = linspace ( a, b, n+1 )
%

u1 = zeros ( 2, n+1 );
u2 = zeros ( 2, n+1 );
x = linspace ( a, b, n+1 );
h = (b-a)/n;

if ( alpha(2) == 0 ) 
   ivp3 = 0;
   u1(:,1) = [alpha(3)/alpha(1); 0];
   u2(:,1) = [0;1];
elseif ( alpha(1) == 0 )
   ivp3 = 0;
   u1(:,1) = [0; alpha(3)/alpha(2)];
   u2(:,1) = [1;0];
else
   ivp3 = 1;
   u3 = zeros ( 2, n+1 );
   u1(:,1) = [0;0];
   u2(:,1) = [1;0];
   u3(:,1) = [0;1];
end;

for i = 1 : n

    [p q r] = feval ( coeff, x(i) );
	k11 = h * [ 0 1; q p ] * u1(:,i) + h * [0; r];
    k12 = h * [ 0 1; q p ] * u2(:,i);
	if ( ivp3 ) k13 = h * [ 0 1; q p ] * u3(:,i); end;
	
	[p q r] = feval ( coeff, x(i) + h/2 );
	k21 = h * [ 0 1; q p ] * ( u1(:,i) + 0.5*k11 ) + h * [0; r];
    k22 = h * [ 0 1; q p ] * ( u2(:,i) + 0.5*k12 );
	if ( ivp3 ) k23 = h * [ 0 1; q p ] * ( u3(:,i) + 0.5*k13 ); end;
	
	k31 = h * [ 0 1; q p ] * ( u1(:,i) + 0.5*k21 ) + h * [0; r];
    k32 = h * [ 0 1; q p ] * ( u2(:,i) + 0.5*k22 );
	if ( ivp3 ) k33 = h * [ 0 1; q p ] * ( u3(:,i) + 0.5*k23 ); end;
	
	[p q r] = feval ( coeff, x(i) + h );
	k41 = h * [ 0 1; q p ] * ( u1(:,i) + k31 ) + h * [0; r];
    k42 = h * [ 0 1; q p ] * ( u2(:,i) + k32 );
	if ( ivp3 ) k43 = h * [ 0 1; q p ] * ( u3(:,i) + k33 ); end;
	
    u1(:,i+1) = u1(:,i) + ( k11 + 2*k21 + 2*k31 + k41 ) / 6;
	u2(:,i+1) = u2(:,i) + ( k12 + 2*k22 + 2*k32 + k42 ) / 6;
	if ( ivp3 ) u3(:,i+1) = u3(:,i) + ( k13 + 2*k23 + 2*k33 + k43 ) / 6; end
	
end;

if ( ~ivp3 )
   if ( beta(2) == 0 )
      c = ( beta(3)/beta(1) - u1(1,n+1) ) / u2(1,n+1);
   elseif ( beta(1) == 0 )
      c = ( beta(3)/beta(2) - u1(2,n+1) ) / u2(2,n+1);
   else
      c = ( beta(3) - beta(1) * u1(1,n+1) - beta(2) * u1(2,n+1) ) / ( beta(1) * u2(1,n+1) + beta(2) * u2(2,n+1) );
   end;
   w = u1 + c * u2;
else
   if ( beta(2) == 0 )
      denom = alpha(1) * u3(1,n+1) - alpha(2) * u2(1,n+1);
	  rhs = beta(3) / beta(1) - u1(1,n+1);
	  c1 = ( alpha(3) * u3(1,n+1) - alpha(2) * rhs ) / denom;
	  c2 = ( alpha(1) * rhs - alpha(3) * u2(1,n+1) ) / denom;
   elseif ( beta(1) == 0 )
      denom = alpha(1) * u3(2,n+1) - alpha(2) * u2(2,n+1);
	  rhs = beta(3) / beta(2) - u1(2,n+1);
	  c1 = ( alpha(3) * u3(2,n+1) - alpha(2) * rhs ) / denom;
	  c2 = ( alpha(1) * rhs - alpha(3) * u2(2,n+1) ) / denom;
   else
      a1 = beta(1) * u2(1,n+1) + beta(2) * u2(2,n+1);
	  a2 = beta(1) * u3(1,n+1) + beta(2) * u3(2,n+1);
	  rhs = beta(3) - beta(1) * u1(1,n+1) - beta(2) * u1(2,n+1);
	  denom = alpha(1) * a2 - alpha(2) * a1;
	  c1 = ( alpha(3) * a2 - alpha(2) * rhs ) / denom;
	  c2 = ( alpha(1) * rhs - alpha(3) * a1 ) / denom;
   end;
   w = u1 + c1 * u2 + c2 * u3;
end;

	
		
