% Some test values
x=1; y=1;
% Peaks function
% z = f(x,y)= 
z=3*(1-x).^2.*exp(-(x.^2) - (y+1).^2) - ...
    10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) ...
    - 1/3*exp(-(x+1).^2 - y.^2)
% Derivatives of the peaks function
% First derivative with respect to x
dzdx = -6*(1-x)*exp(-x^2-(y+1)^2) - ...
         6*(1-x)^2*x*exp(-x^2-(y+1)^2) - ...
         10*(1/5-3*x^2)*exp(-x^2-y^2) + ...
         20*(1/5*x-x^3-y^5)*x*exp(-x^2-y^2) - ...
         1/3*(-2*x-2)*exp(-(x+1)^2-y^2)
% First derivative with respect to y
dzdy = 3*(1-x)^2*(-2*y-2)*exp(-x^2-(y+1)^2) + ...
        50*y^4*exp(-x^2-y^2) + ...
        20*(1/5*x-x^3-y^5)*y*exp(-x^2-y^2) + ...
        2/3*y*exp(-(x+1)^2-y^2)
% 2nd derivative 
ddzdxdx = 36*x*exp(-x^2-(y+1)^2) - 18*x^2*exp(-x^2-(y+1)^2) - ...
              24*x^3*exp(-x^2-(y+1)^2) + 12*x^4*exp(-x^2-(y+1)^2) + ...
              72*x*exp(-x^2-y^2) - 148*x^3*exp(-x^2-y^2) - ...
              20*y^5*exp(-x^2-y^2) + 40*x^5*exp(-x^2-y^2) + ...
              40*x^2*exp(-x^2-y^2)*y^5 -2/3*exp(-(x+1)^2-y^2) - ...
              4/3*exp(-(x+1)^2-y^2)*x^2 -8/3*exp(-(x+1)^2-y^2)*x
ddzdydy = -6*(1-x)^2*exp(-x^2-(y+1)^2) + 3*(1-x)^2*(-2*y-2)^2*exp(-x^2-(y+1)^2) + ...
          200*y^3*exp(-x^2-y^2)-200*y^5*exp(-x^2-y^2) + ...
          20*(1/5*x-x^3-y^5)*exp(-x^2-y^2) - 40*(1/5*x-x^3-y^5)*y^2*exp(-x^2-y^2) + ...
          2/3*exp(-(x+1)^2-y^2)-4/3*y^2*exp(-(x+1)^2-y^2)    
% ddzdxdy = 12*(y+1)*(1-x)*exp(-x^2-(y+1)^2) + ...
%          12*(y+1)*(1-x)^2*x*exp(-x^2-(y+1)^2) + ...
%          20*y*(1/5-3*x^2)*exp(-x^2-y^2) - ...
%          40*y*(1/5*x-x^3-y^5)*x*exp(-x^2-y^2) + ...
%          100*y^4*x*exp(-x^2-y^2) + ...
%          2/3*y*(-2*x-2)*exp(-(x+1)^2-y^2)
% ddzdydx =-6*x*(1-x)^2*(-2*y-2)*exp(-x^2-(y+1)^2) - ...
%         6*(1-x)*(-2*y-2)*exp(-x^2-(y+1)^2) - ...
%         100*x*y^4*exp(-x^2-y^2) + ...
%         4*y*exp(-x^2-y^2) - ...
%         60*x^2*y*exp(-x^2-y^2) - ...
%         40*x*(1/5*x-x^3-y^5)*y*exp(-x^2-y^2) - ...
%         4/3*(x+1)*y*exp(-(x+1)^2-y^2)
