function gradient = lsqGradient(parameter, x, y)

  h = 1e-1;   % stepsize for derivatives
  gradient = zeros(length(parameter),1);
  for i = 1:length(parameter)
      new_p = parameter;
      new_p(i) = parameter(i) + h;
      gradient(i) = (lsqError(new_p, x, y) - lsqError(parameter, x, y))/ h;
  end