function [ A, b ] = plasma_matrix ( n )

%*****************************************************************************80
%
%% plasma_matrix() computes jacobian and residual for a plasma problem.
%
%  Licensing:
%
%    This code is distributed under the GnU LGPL license.
%
%  Modified:
%
%    10 March 2020
%
%  Author:
%
%    James Cheung;
%    Modifications by John Burkardt.
%
%  Input:
%
%    integer N: the number of nodes in the X and Y directions.
%
%  Output:
%
%    real sparse A(n*n,n*n): the Jacobian matrix.
%
%    real b(n*n,1): the residual vector.
%
%  Local:
%
%    real XLEFT, XRIGHT: the left and right endpoints of the X and Y intervals.
%
%    real H: the spacing between nodes.
%
%    real RHO(n,n): the density. 
%
%    integer NUMNODES: the number of nodes.
%
%    real X(n): the X coordinates of the grid.
%
%    real Y(n): the Y coordinates of the grid.
%
  xleft = -10.0;
  xright = +10.0;

  if ( nargin < 1 )
    n = 101;
  end

  x = ( linspace ( xleft, xright, n ) )';
  y = x;
  h = x(2) - x(1);
%
%  Set up the plasma density.
%
  numnodes = n * n;
  rho = zeros ( numnodes, 1 );
  phi = zeros ( numnodes, 1 );
%
%  Note that, initially, terms like
%  - h^2 * exp ( phi(*) ) + h^2 rho(*) are zero.
%  In a fuller version of the code, these terms drift apart.
%
  nr = 0.25;
  k = 0;

  for i = 1 : n
    for j = 1 : n

      k = k + 1;

      if ( abs ( x(i) ) < 5.0 && abs ( y(j) ) < 5.0 )
        rho(k) = 1.0;
      else
        rho(k) = nr;
      end

      phi(k) = log ( rho(k) );

    end
  end
%
%  Allocate space.
%
  A = spalloc ( numnodes, numnodes, 5*numnodes );
  b = zeros ( numnodes, 1 );
%
%  Bottom
%
  b(1) = ...
    - 4.0 * phi(1) ...
    + 2.0 * phi(2) ...
    + 2.0 * phi(n+1) ...
    - h^2 * exp ( phi(1) ) ...
    + h^2 * rho(1);

  A(1,1) = - 4.0 - h^2 * exp ( phi(1) );
  A(1,2) = 2.0;
  A(1,n+1) = 2.0;

  for i = 2 : n - 1

    b(i) = ...
              phi(i-1) ...
      - 4.0 * phi(i) ...
      +       phi(i+1) ...
      + 2.0 * phi(i+n) ...
      - h^2 * exp ( phi(i) ) ...
      + h^2 * rho(i);

    A(i,i-1) = 1.0;
    A(i,i) = - 4.0 - h^2 * exp ( phi(i) );
    A(i,i+1) = 1.0;
    A(i,n+i) = 2.0;

  end

  b(n) = ...
      2.0 * phi(n-1) ...
    - 4.0 * phi(n) ...
    + 2.0 * phi(2*n) ...
    - h^2 * exp ( phi(n) ) ...
    + h^2 * rho(n); 

  A(n,n-1) = 2.0;
  A(n,n) = - 4.0 - h^2 * exp ( phi(n) );
  A(n,n+n) = 2.0;
%
%  Middle
%
  for i = 2 : n - 1

    z = ( i - 1 ) * n;

    b(z+1) = ...
              phi(z+1-n) ...
      - 4.0 * phi(z+1) ...
      + 2.0 * phi(z+2) ...
      +       phi(z+1+n) ...
      - h^2 * exp ( phi(z+1) ) ...
      + h^2 * rho(z+1);

    A(z+1,z-n+1) = 1.0;
    A(z+1,z+1) = - 4.0 - h^2 * exp ( phi(z+1) );
    A(z+1,z+2) = 2.0;
    A(z+1,z+n+1) = 1.0;

    for j = 2 : n - 1

      b(z+j) = ...
                phi(z+j-1) ...
        +       phi(z+j-n) ...
        - 4.0 * phi(z+j) ...
        +       phi(z+j+1) ...
        +       phi(z+n+j) ...
        - h^2 * exp ( phi(z+j) ) ...
        + h^2 * rho(z+j);

      A(z+j,z-n+j) = 1.0;
      A(z+j,z+j-1) = 1.0;
      A(z+j,z+j) = - 4.0 - h^2 * exp ( phi(z+j) );
      A(z+j,z+j+1) = 1.0;
      A(z+j,z+n+j) = 1.0;

    end

    b(z+n) = ...
        2.0 * phi(z+n-1) ...
      +       phi(z) ...
      - 4.0 * phi(z+n) ...
      +       phi(z+2*n) ...
      - h^2 * exp ( phi(z+n) ) ...
      + h^2 * rho(z+n);

    A(z+n,z) = 1.0;
    A(z+n,z+n-1) = 2.0;
    A(z+n,z+n) = - 4.0 - h^2 * exp ( phi(z+n) );
    A(z+n,z+2*n) = 1.0;

  end
%
%  Top
%
  z = ( n - 1 ) * n;

  b(z+1) = ...
      2.0 * phi(z-n+1) ...
    - 4.0 * phi(z+1) ...
    + 2.0 * phi(z+2) ...
    - h^2 * exp ( phi(z+1) ) ...
    + h^2 * rho(z+1);

  A(z+1,z-n+1) = 2.0;
  A(z+1,z+1) = - 4.0 - h^2 * exp ( phi(z+1) );
  A(z+1,z+2) = 2.0;

  for i = 2 : n - 1

    b(z+i) = ...
              phi(z+i-1) ...
      + 2.0 * phi(z+i-n) ...
      - 4.0 * phi(z+i) ...
      + phi(z+i+1) ...
      - h^2 * exp ( phi(z+i) ) ...
      + h^2 * rho(z+i);

    A(z+i,z-n+i) = 2.0;
    A(z+i,z+i-1) = 1.0;
    A(z+i,z+i) = - 4.0 - h^2 * exp ( phi(z+i) );
    A(z+i,z+i+1) = 1.0;

  end

  b(z+n) = ...
      2.0 * phi(z+n-1) ...
    + 2.0 * phi(z) ... 
    - 4.0 * phi(z+n) ...
    - h^2 * exp ( phi(z+n) ) ...
    + h^2 * rho(z+n);

  A(z+n,z) = 2.0;
  A(z+n,z+n-1) = 2.0;
  A(z+n,z+n) = - 4.0 - h^2 * exp ( phi(z+n) );

  return
end
 
