function fv=parab1(v0)
a=[1;1;1];
v0=1.15*(rand(3,1)-a/2) % starting values are generated randomly
     h = .5;            % step size
  v=v0;
ff=f(v0);
step=0;
while(norm(ff)>1e-12)
     jj=[f(v+[h;0;0])-f(v-[h;0;0]),f(v+[0;h;0])-f(v-[0;h;0]),f(v+[0;0;h])-f(v-[0;0;h])]; % jacobian
     delta=-jj\ff;
     v=v+delta;  % quasi-Newton
     ff=f(v);
 step=step+1;
end

 step, v,
function ff=f(v)
ff=[v'*v-1; (v(1)-1/2)^2/3+v(2)^2+3*v(3)^2+v(1)*v(3)-1; v(1)^2/3+3*(v(2)+1/4)^2+2*v(3)^2-v(2)*v(3)+3/2*v(1)*v(2)-1];

step = 6
v = 0.4089
   -0.8899
    0.2022
step = 9
v = -0.8713
    -0.4759
    -0.1193
step = 6
v = 0.7059
   -0.1135
   -0.6992
step = 11
v =-0.8474
   -0.0283
    0.5302
step = 8
v =-0.9210
    0.3618
   -0.1447
step = 7
v = 0.9004
    0.0658
    0.4301
step = 8
v = 0.2492
   -0.9435
   -0.2186
step = 6
v =-0.8376
    0.1676
    0.5199
