I am replicating the following equation 7 from *Daria Nesterovich Anderson et al 2018 J. Neural Eng. 15 026005*:

When I implement my code, my optimization is always the same no matter what I use as my target region, omega. Any thoughts on where I’ve gone wrong? Thank you.

My code is below:

```
size_omega = sum(sum(sum(omega)));
Integral_omega = zeros(1,8);
for i=1:8
if (i==1)
[contact1] = voltagebox(contact1);
[all_derivatives] = Part2ndDeriv(contact1);
dy2 = all_derivatives{2,1};
dx2 = all_derivatives{2,2};
dxy = all_derivatives{2,4};
dyx = all_derivatives{2,7};
dz2 = all_derivatives{2,3};
dxz = all_derivatives{2,5};
dzx = all_derivatives{2,8};
dyz = all_derivatives{2,6};
dzy = all_derivatives{2,9};
elseif (i==2)
[contact2] = voltagebox(contact2);
[all_derivatives] = Part2ndDeriv(contact2);
dy2 = all_derivatives{2,1};
dx2 = all_derivatives{2,2};
dxy = all_derivatives{2,4};
dyx = all_derivatives{2,7};
dz2 = all_derivatives{2,3};
dxz = all_derivatives{2,5};
dzx = all_derivatives{2,8};
dyz = all_derivatives{2,6};
dzy = all_derivatives{2,9};
elseif (i==3)
[contact3] = voltagebox(contact3);
[all_derivatives] = Part2ndDeriv(contact3);
dy2 = all_derivatives{2,1};
dx2 = all_derivatives{2,2};
dxy = all_derivatives{2,4};
dyx = all_derivatives{2,7};
dz2 = all_derivatives{2,3};
dxz = all_derivatives{2,5};
dzx = all_derivatives{2,8};
dyz = all_derivatives{2,6};
dzy = all_derivatives{2,9};
elseif (i==4)
[contact4] = voltagebox(contact4);
[all_derivatives] = Part2ndDeriv(contact4);
dy2 = all_derivatives{2,1};
dx2 = all_derivatives{2,2};
dxy = all_derivatives{2,4};
dyx = all_derivatives{2,7};
dz2 = all_derivatives{2,3};
dxz = all_derivatives{2,5};
dzx = all_derivatives{2,8};
dyz = all_derivatives{2,6};
dzy = all_derivatives{2,9};
elseif (i==5)
[contact5] = voltagebox(contact5);
[all_derivatives] = Part2ndDeriv(contact5);
dy2 = all_derivatives{2,1};
dx2 = all_derivatives{2,2};
dxy = all_derivatives{2,4};
dyx = all_derivatives{2,7};
dz2 = all_derivatives{2,3};
dxz = all_derivatives{2,5};
dzx = all_derivatives{2,8};
dyz = all_derivatives{2,6};
dzy = all_derivatives{2,9};
elseif (i==6)
[contact6] = voltagebox(contact6);
[all_derivatives] = Part2ndDeriv(contact6);
dy2 = all_derivatives{2,1};
dx2 = all_derivatives{2,2};
dxy = all_derivatives{2,4};
dyx = all_derivatives{2,7};
dz2 = all_derivatives{2,3};
dxz = all_derivatives{2,5};
dzx = all_derivatives{2,8};
dyz = all_derivatives{2,6};
dzy = all_derivatives{2,9};
elseif (i==7)
[contact7] = voltagebox(contact7);
[all_derivatives] = Part2ndDeriv(contact7);
dy2 = all_derivatives{2,1};
dx2 = all_derivatives{2,2};
dxy = all_derivatives{2,4};
dyx = all_derivatives{2,7};
dz2 = all_derivatives{2,3};
dxz = all_derivatives{2,5};
dzx = all_derivatives{2,8};
dyz = all_derivatives{2,6};
dzy = all_derivatives{2,9};
elseif (i==8)
[contact8] = voltagebox(contact8);
[all_derivatives] = Part2ndDeriv(contact8);
dy2 = all_derivatives{2,1};
dx2 = all_derivatives{2,2};
dxy = all_derivatives{2,4};
dyx = all_derivatives{2,7};
dz2 = all_derivatives{2,3};
dxz = all_derivatives{2,5};
dzx = all_derivatives{2,8};
dyz = all_derivatives{2,6};
dzy = all_derivatives{2,9};
end
for x=1:199
for y=1:199
for z=1:199
if (omega(x,y,z)==1) %if (x,y,z) is a point in the ROI omega:
% Find the hessian at (x,y,z):
hess = [dx2(x,y,z),dxy(x,y,z),dxz(x,y,z); dyx(x,y,z),dy2(x,y,z),dyz(x,y,z); dzx(x,y,z),dzy(x,y,z),dz2(x,y,z)];
a=(isnan(hess) | isinf(hess));
if (sum(sum(a)>0))
[m n] = find(a==1);
aa = length(m);
for b=1:aa
hess(m(b),n(b)) = 0;
end
end
% Determine the eigenvalues of the hessian:
[~, val] = eig(hess);
% Sum the eigenvalues:
lambda_sum = val(1)+val(2)+val(3);
Integral_omega(i) = lambda_sum + Integral_omega(i);
else
Integral_omega(i) = Integral_omega(i) + 0;
end
end
end
end
end
cvx_begin
variables c(i)
maximize sum(c(i)/size_omega * Integral_omega)
subject to
c<=4
cvx_end
```