chandra.m

linear equation routines in matlab

function varargout = chandra(x); 
%
%  Compute the residual of the discretized 
%  Chandrasekhar's H-equation
%
global N c


mmu = zeros(N,1);
for i = 1:N
    mmu(i) = (i-0.5)/N;
end
  
  

% return the function value
F = x;
for i = 1:N
     tmp = sum( mmu(i)*x ./ (mmu(i) + mmu) );
     F(i) = F(i) - 1/(1 - c*tmp/(2*N) );
end
varargout{1} = F;

if( nargout > 1 )
  % return the Jacobian as the second argument
  Jac = eye(N);
  for i = 1:N
     tmp  = sum( mmu(i)*x ./ (mmu(i) + mmu) );
     for j = 1:N
         Jac(i,j) = Jac(i,j) - ...
                    (c/(2*N))*mmu(i)/(mmu(i) + mmu(j))/(1 - c*tmp/(2*N) )^2;
     end
  end
  varargout{2} = Jac;
end