function ill_bvp_bvp4c ( epsilon )

%*****************************************************************************80
%
%% ill_bvp_bvp4c() calls bvp4c() to solve the ill conditioned BVP.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    21 June 2022
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real EPSILON, the value of the parameter.
%

%
%  Compute an initial solution estimate.
%
  a = -1.0;
  b = +1.0;
  x_init = linspace ( a, b, 5 );
  y_init = [ 2.0, 0.0 ];
  solinit = bvpinit ( x_init, y_init );
%
%  Request a solution from bvp4c():
%
  sol = bvp4c ( @ill_bvp_ode, @ill_bvp_bc, solinit );
%
%  Evaluate the solution at 101 plot points.
%
  x = linspace ( a, b, 101 );
  y = deval ( sol, x );
%
%  Display a plot of the solution.
%
  clf ( );
  plot ( x, y(1,:), 'r-', 'Linewidth', 2 );
  xlabel ( '<--- X --->' );
  ylabel ( '<--- Y --->' );
  title ( sprintf ( 'ill conditioned BVP, \\epsilon = %g\n', epsilon ) );
  grid ( 'on' );
  pbaspect ( [3,1,1] );
% axis ( 'equal' );
  filename = sprintf ( 'ill_bvp_%f.png', epsilon );
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s"\n', filename );

  return
end
function dydx = ill_bvp_ode ( x, y )

%*****************************************************************************80
%
%% ill_bvp_ode() evaluates the right hand side of the ill conditioned BVP.
%
%  Discussion:
%
%    We assume that the differential equation has been rewritten as a
%    system of first order equations of the form
%
%      dydx = f(x,y)
%
%    The second order system:
%
%      eps y'' - xy' + y = 0
%
%    is rewritten as
%
%      u' = v
%      v' = ( xv - u ) / eps  
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    21 June 2022
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real X, the point at which the ODE is to be evaluated.
%
%    real Y(M), the value of the solution at X.
%
%  Output:
%
%    real DYDX(M), the value of the right hand side given X and Y.
%
  global epsilon

  dydx(1) = y(2);
  dydx(2) = ( x * y(2) - y(1) ) / epsilon;

  return
end
function bc = ill_bvp_bc ( ya, yb )

%*****************************************************************************80
%
%% ill_bvp_bc() evaluates the boundary conditions for the ill conditioned bvp.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    21 June 2022
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    real YA(M), YB(M), the solution value at the left and right endpoints.
%
%  Output:
%
%    real BC(2), the value of the boundary conditions.
%
  bc(1) = ya(1) - 2.0;
  bc(2) = yb(1) - 1.0;

  return
end

