function d=grad(fname,x,param,delta,ind);
%GRAD    Compute gradients
%        d=grad(func,x,param,delta,ind) takes as input the name of a 
%        function (func) which is a string and the vectors x and param
%        which will be the inputs for func.m (e.g. y=func(x,param)).  
%        func.m has as output some matrix Y (=func(x,param)). grad 
%        computes d=d(vec(Y))/dx with deviations from delta. The (i,j)
%        element is the derivative of the ith function with respect
%        to the jth element of x.  delta is a vector with the same 
%        length as x.  d=grad(func,x,param) uses delta=abs(x)*1e-4 
%        and ind=[1:length(x)]'.

%        Ellen R. McGrattan,  2-17-90
%        revised, 2-13-93

n=length(x);
if nargin < 5;
  ind=[1:n]';
end;
m=length(ind);
if nargin < 4;
  delta=abs(x)*1e-4; 
end;
if n~=length(delta);
  del=zeros(x);
  del(ind)=delta;
  delta=del;
end;
delta=diag(max(delta,1e-8));
for j=1:m;
  k=ind(j);
  eval(['d(:,j)=vec( ',fname,'(x+delta(:,k),param)-',fname, ..
        '(x-delta(:,k),param) )/(2*delta(k,k)) ;']);
end;

