octave:8> D1=rand(4)
D1 =

   0.68272   0.98911   0.82106   0.96389
   0.67374   0.79376   0.17256   0.14950
   0.59830   0.53012   0.77591   0.86876
   0.84454   0.25435   0.25570   0.99262

octave:9> F=@(u) (D1*u).^2;
octave:10> JF=@(u) 2*diag(D1*u)*D1;
octave:11> epsilon=1e-4;
octave:12> u=rand(4,1)
u =

   0.92854
   0.53687
   0.76574
   0.23354

octave:13> v=rand(4,1)
v =

   0.79954
   0.53031
   0.21695
   0.56006

octave:14> (F(u+eps*v)-F(u))/eps
ans =

   8
   5
   4
   3

octave:15> (F(u+epsilon*v)-F(u))/epsilon
ans =

   7.2210
   2.6346
   4.6314
   3.8337

octave:16> JF(u)*v
ans =

   7.2206
   2.6344
   4.6312
   3.8335

octave:17> JF=@(u) D1^2*u.^u;
octave:18> JF(u)*v
error: operator *: nonconformant arguments (op1 is 4x1, op2 is 4x1)
octave:18> JF=@(u) D1^2*u.^2;
octave:19> JF(u)*v
error: operator *: nonconformant arguments (op1 is 4x1, op2 is 4x1)
octave:19> JF=@(u) D1^2*(u.^2);
octave:20> JF(u)*v
error: operator *: nonconformant arguments (op1 is 4x1, op2 is 4x1)
octave:20> JF=@(u) (D1^2)*(u.^2);
octave:21> JF(u)*v
error: operator *: nonconformant arguments (op1 is 4x1, op2 is 4x1)
octave:21> D1
D1 =

   0.68272   0.98911   0.82106   0.96389
   0.67374   0.79376   0.17256   0.14950
   0.59830   0.53012   0.77591   0.86876
   0.84454   0.25435   0.25570   0.99262

octave:22> D2^2
error: 'D2' undefined near line 1 column 1
octave:22> D1^2
ans =

   2.43780   2.14082   1.61478   2.47602
   1.22425   1.42594   0.86226   1.06637
   1.96356   1.64485   1.40690   2.19237
   1.73925   1.42525   1.18953   2.05951

octave:23> F=@(u) = D1*(u.^2);
parse error:

  syntax error

>>> F=@(u) = D1*(u.^2);
           ^

octave:23> F=@(u) D1*(u.^2);
octave:24> JF=@(u) 2*D1*diag(u);
octave:25> (F(u+epsilon*v)-F(u))/epsilon
ans =

   2.1020
   1.5488
   1.6754
   1.7435

octave:26> JF=@(u) 2*D1*diag(u);
octave:27> JF(u)*v
ans =

   2.1019
   1.5488
   1.6753
   1.7434

octave:28> JF=@(u) 2*diag(u)*D1;
octave:29> JF(u)*v
ans =

   3.32113
   1.16047
   2.16609
   0.66397

octave:30> F=@(u) 2*u.*(D1*u);
octave:31> JF=@(u) 2*diag(D1*u)+2*diag(u)*D1;
octave:32> (F(u+epsilon*v)-F(u))/epsilon
ans =

   6.5496
   2.4532
   2.8765
   2.1745

octave:33> JF(u)*v
ans =

   6.5493
   2.4531
   2.8765
   2.1743

octave:34> diary off