function d=grad(fname,x,delta,ind,par);
%GRAD    Compute gradients
%        d=grad(func,x,delta,ind) takes as input the name of a function
%        (func) which is a string and a vector x which will be the
%        input for func.m.  func.m has as output some matrix Y 
%        (=func(x)). grad computes d=d(vec(Y))/dx with deviations 
%        from delta.  delta is a vector with the same length as x.
%        d=grad(func,x) uses delta=abs(x)*1e-4 and ind=[1:length(x)]'.

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

if nargin < 5;
  par = 0;
end;
n=length(x);
if nargin < 4;
  ind=[1:n]';
end;
m=length(ind);
if nargin < 3;
  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),par)-',fname, ...
        '(x-delta(:,k),par) )/(2*delta(k,k)) ;']);
end;

