  function[alpha,beta,chisq] = mrqcof(x,y,sig,ndata,ma,a);

  alpha=zeros(ma,ma);
  beta=zeros(1,ma);
  g=zeros(1,ndata);
  g=1.0./sig.^2;

  chisq=0;

  [ymod,dyda] = funcs(x,a,ma,ndata);

  dy=y-ymod;
  for k=1:ndata
    for i=1:ma
      for j=1:ma
        alpha(i,j)=alpha(i,j)+ g(k)*dyda(i,k)*dyda(j,k);
      end
      beta(i)=beta(i)+g(k)*dy(k)*dyda(i,k);
    end
    chisq=chisq+g(k)*dy(k)^2;
  end

