function stab ( )

%*****************************************************************************80
%
%% stab() examines mean-square and asymptotic stability.
%
%  Discussion:
%
%    The function tests the mean-square and asymptotic stability
%    of the Euler-Maruyama method applied to a stochastic differential
%    equation (SDE).
%
%    The SDE is
%
%      dX = lambda*X dt + mu*X dW,
%      X(0) = Xzero,
%
%    where 
%
%      lambda is a constant,
%      mu is a constant,
%      Xzero = 1.
%
%    We will examine different values of lambda and mu.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    04 October 2006
%
%  Author:
%
%    Desmond Higham
%
%  Reference:
%
%    Desmond Higham,
%    An Algorithmic Introduction to Numerical Simulation of
%    Stochastic Differential Equations,
%    SIAM Review,
%    Volume 43, Number 3, September 2001, pages 525-546.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'stab() checks asymptotic and mean-square stability.\n' );

% rng ( 'shuffle' );
%
%  Set problem parameters.
%
  tmax = 20.0;
  m = 50000;
  xzero = 1;          
%
%  Set line types for the plot.
%
  ltype = { 'b-', 'r--', 'm-.' };
%
%  Mean square.
%
  subplot ( 211 );

  lambda = -3.0;
  mu = sqrt ( 3.0 );
%
%  XMS is the mean square estimate of M paths.
%
  hold ( 'on' );
  for k = 1 : 3

    dt = 2^(1-k);                      
    n = tmax / dt;
    xms = zeros(1,n);
    xtemp = xzero * ones(m,1);

    for j = 1 : n
      winc = sqrt ( dt ) * randn(m,1);  
      xtemp = xtemp + dt * lambda * xtemp + mu * xtemp .* winc;
      xms(j) = mean ( xtemp.^2 );
    end

    semilogy ( [0:dt:tmax], [xzero,xms], ltype{k}, 'Linewidth', 2 ), 

  end

  legend ( '\Delta t = 1','\Delta t = 1/2','\Delta t = 1/4' );
  title ( 'Mean-Square: \lambda = -3, \mu = \surd 3', 'FontSize', 16 );
  ylabel ( 'E[X^2]','FontSize', 12 );
  axis ( [0,tmax,1e-20,1e+20] );
  hold ( 'off' );
%
%  Asymptotic, a single path.
%
  subplot ( 212 );
  tmax = 500.0;
%
%  Problem parameters.
%
  lambda = 0.5;
  mu = sqrt ( 6.0 );

  hold ( 'on' );

  for k = 1 : 3

    dt = 2^(1-k);                     
    n = tmax / dt;
    xemabs = zeros(1,n);
    xtemp = xzero;

    for j = 1 : n
      winc = sqrt ( dt ) * randn;      
      xtemp = xtemp + dt * lambda * xtemp + mu * xtemp * winc;
      xemabs(j) = abs ( xtemp );
    end

    semilogy ( [0:dt:tmax], [xzero,xemabs], ltype{k}, 'Linewidth', 2 );

  end

  legend ( '\Delta t = 1','\Delta t = 1/2','\Delta t = 1/4' );
  title ( 'Single Path: \lambda = 1/2, \mu = \surd 6','FontSize',16 );
  ylabel ( '|X|','FontSize',12 );
  axis ( [0,tmax,1e-50,1e+100] );
  hold ( 'off' );

  filename = 'stab.png';
  print ( '-dpng', filename );
  fprintf ( 1, '  Graphics saved as "%s".\n', filename );

  return
end
