function solinit = bvpinit(x,v,parameters)
%BVPINIT  Form the initial guess for BVP4C.
%   SOLINIT = BVPINIT(X,YINIT) forms the initial guess for BVP4C in common
%   circumstances. The BVP is to be solved on [a,b]. The vector X specifies
%   a and b as X(1) = a and X(end) = b. It is also a guess for an appropriate
%   mesh. BVP4C will adapt this mesh to the solution, so often a guess like
%   X = linspace(a,b,10) will suffice, but in difficult cases, mesh points must
%   be placed where the solution changes rapidly. 
%
%   The entries of X must be ordered and distinct, so if a < b, then
%   X(1) < X(2) < ... < X(end), and similarly for a > b. 
%
%   YINIT provides a guess for the solution. It must be possible to evaluate the
%   differential equations and boundary conditions for this guess. YINIT can be
%   either a vector or a function:
%
%   vector:  YINIT(i) is a constant guess for the i-th component Y(i,:) of the
%            solution at all the mesh points in X.
%
%   function:  YINIT is a function of a scalar x. For example, use 
%              solinit = bvpinit(x,@yfun), if for any x in [a,b], yfun(x) 
%              returns a guess for the solution y(x).
%                       
%   SOLINIT = BVPINIT(X,YINIT,PARAMETERS) indicates that the BVP involves unknown
%   parameters. A guess must be provided for all parameters in the vector
%   PARAMETERS. 
%
%   See also BVPGET, BVPSET, BVPVAL, BVP4C, @.

%   Jacek Kierzenka and Lawrence F. Shampine
%   Copyright 1984-2000 The MathWorks, Inc. 
%   $Revision: 1.6 $  $Date: 2000/06/02 00:11:27 $

N = length(x);
if x(1) < x(N)
  if any(diff(x) <= 0)
    error('The entries of x must satisfy a = x(1) < x(2) < ... < x(end) = b.')
  end
else
  if any(diff(x) >= 0)
    error('The entries of x must satisfy a = x(1) > x(2) > ... > x(end) = b.')
  end
end

if isnumeric(v) 
  w = v;
else
  w = feval(v,x(1));
end
[m,n] = size(w);
if m == 1
  L = n;
elseif n == 1
  L = m;
else
  error('The guess for the solution must return a vector.')
end

yinit = zeros(L,N);
if isnumeric(v)
  yinit = repmat(v(:),1,N);
else 
  yinit(:,1) = w(:);
  for i = 2:N
    w = feval(v,x(i));
    yinit(:,i) = w(:);
  end
end

solinit.x = x(:)';
solinit.y = yinit;
if nargin > 2 
  solinit.parameters = parameters;
end
