gradient descent algorithm - index out of bounds

I have written a simple code for the gradient descent algorithm but i am getting an error while accessing grd(2) . It says " Attempted to access grd(2); index out of bounds because numel(grd)=1." . I am not able to find out the bug. If anyone expert in it ,. can you please tell me my where i am doing wrong.

    function [xopt,fopt,niter,gnorm,dx] = mystochastic(varargin)

      if nargin==0
          % define starting point
          x0 = [5 5]';
      elseif nargin==1
          % if a single input argument is provided, it is a user-defined starting
          % point.
          x0 = varargin{1};
      else
          error('Incorrect number of input arguments.')
      end

      % termination tolerance
      tol = 1e-6;

      % maximum number of allowed iterations
      maxiter = 1000;

      % minimum allowed perturbation
      dxmin = 1e-6;

      % step size ( 0.33 causes instability, 0.2 quite accurate)
      alpha = 0.1;

      % initialize gradient norm, optimization vector, iteration counter, perturbation
      gnorm = inf; x = x0; niter = 0; dx = inf;

      % define the objective function:
      f = @(x1,x2) 1.7*exp(-(((x1-3).^2)/10 + ((x2-3).^2)/10)) + exp(-(((x1+5).^2)/8 + ((x2+5).^2)/8))+2*exp(-((x1.^2)/4 + ((x2.^2)/5)))+ 1.5*exp(-(((x1-4).^2)/18 + ((x2+4).^2)/16)) + 1.2 *exp(-(((x1+4).^2)/18 + ((x2-4).^2)/16));

      % plot objective function contours for visualization:
      figure(1); clf; ezcontour(f,[-10 10 -10 10]); axis equal; hold on

      % gradient descent algorithm:
      while and(gnorm>=tol, and(niter <= maxiter, dx >= dxmin))
          % calculate gradient:
          r = randi([1 2]);
          g = grad(x,r);
          gnorm = norm(g);
          % take step:
          xnew = x - alpha*g;
          % check step
          if ~isfinite(xnew)
              display(['Number of iterations: ' num2str(niter)])
              error('x is inf or NaN')
          end
          % plot current point
          plot([x(1) xnew(1)],[x(2) xnew(2)],'ko-')
          refresh
          % update termination metrics
          niter = niter + 1;
          dx = norm(xnew-x);
          x = xnew;

      end
      xopt = x;
      fopt = f2(xopt);
      niter = niter - 1;

      function g = grad(x,r)
      grd = (- x(1).*exp(- x(1).^2/4 - x(2).^2/5) - (4*x(2).*exp(- x(1).^2/4 - x(2).^2/5))/5 - (17*exp(- (x(1) - 3).^2/10 - (x(2) - 3).^2/10)*(x(1)/5 - 3/5))/10 - exp(- (x(1) + 5).^2/8 - (x(2) + 5).^2/8)*(x(1)/4 + 5/4) - (3*exp(- (x(1) - 4).^2/18 - (x(2) + 4).^2/16)*(x(1)/9 - 4/9))/2 - (6*exp(- (x(1) + 4).^2/18 - (x(2) - 4).^2/16)*(x(1)/9 + 4/9))/5 - (17*exp(- (x(1) - 3).^2/10 - (x(2) - 3).^2/10)*(x(2)/5 - 3/5))/10 - exp(- (x(1) + 5).^2/8 - (x(2) + 5).^2/8)*(x(2)/4 + 5/4) - (3*exp(- (x(1) - 4).^2/18 - (x(2) + 4).^2/16)*(x(2)/8 + 1/2))/2 - (6*exp(- (x(1) + 4).^2/18 - (x(2) - 4).^2/16)*(x(2)/8 - 1/2))/5);
      g = zeros(size(x,1),1);
      g(r) = grd(r);

Answer 1

[Walter Roberson]

        grd = (- x(1).*exp(- x(1).^2/4 - x(2).^2/5) - (4*x(2).*exp(- x(1).^2/4 - x(2).^2/5))/5 - (17*exp(- (x(1) - 3).^2/10 - (x(2) - 3).^2/10)*(x(1)/5 - 3/5))/10 - exp(- (x(1) + 5).^2/8 - (x(2) + 5).^2/8)*(x(1)/4 + 5/4) - (3*exp(- (x(1) - 4).^2/18 - (x(2) + 4).^2/16)*(x(1)/9 - 4/9))/2 - (6*exp(- (x(1) + 4).^2/18 - (x(2) - 4).^2/16)*(x(1)/9 + 4/9))/5 - (17*exp(- (x(1) - 3).^2/10 - (x(2) - 3).^2/10)*(x(2)/5 - 3/5))/10 - exp(- (x(1) + 5).^2/8 - (x(2) + 5).^2/8)*(x(2)/4 + 5/4) - (3*exp(- (x(1) - 4).^2/18 - (x(2) + 4).^2/16)*(x(2)/8 + 1/2))/2 - (6*exp(- (x(1) + 4).^2/18 - (x(2) - 4).^2/16)*(x(2)/8 - 1/2))/5);

never uses [] or cat(), and only refers to x with a scalar subscript, and never uses any other variable. Therefore result is always going to be a scalar, and so the result cannot be indexed at location 2.