Нахождения экстремума (исправить ошибки и добавить матрицу Гессе)

/*Данная программа вычисляет гессиан системы функций. */
/* Используется библиотека ublas */
 
#include "ublas/matrix.hpp"
#include "ublas/io.hpp"
#include "ublas/vector.hpp"
#include <iostream>
#include <cmath>
#include <fstream>
 
//using namespace std;
using namespace boost::numeric::ublas;
//using namespace std;
 
typedef  double (*Function)(vector<double>) ;
 
 
double f(vector<double> point)
{
    //return 2*point[0] * point[0] +point[0] * point[0] +  3*point[1] * point[1];
    return point[0] *point[0] * sqrt(abs(point[1])) ;
}
 
 
 
matrix<double> ComputeHessian( Function f, vector<double> args, double step)
{
    //f - указателей на функцию, для которой вычисляется гессиан
    //args - точка, в которой вычисляется гессиан
    //step - шаг метода вычисления
    //Для представления матрицы используется класс matrix из библиотеки ublas
    
    
    int varCount = args.size();
    matrix<double> hessian(varCount, varCount);
 
    //Вспомогательные 
    vector<double> fPlus(varCount);
    vector<double> fMinus(varCount);
    
 
    double curFValue = f(args); //значение функции в текущей точке 
    for (int i = 0; i < varCount; i++)
    {
        
        /* В случае производной по одной переменной берём в качестве аппроксимации вторую центральную производную */
        fPlus[i] = f(args + step * unit_vector<double>(varCount, i) );
        fMinus[i] = f(args - step * unit_vector<double>(varCount, i) );
 
        hessian(i, i) = (fPlus[i] - 2 * curFValue + fMinus[i]) / (step * step);
    }
    
    for (int i = 0; i < varCount; i++)
        //Так как гессиан симметричен, вычисляем только его половину
        for (int j = i + 1; j < varCount; j++)
        {
 
            double x =  f(args +  step * (unit_vector<double>(varCount, i) + unit_vector<double>(varCount, j) ) );
            hessian(i, j) =  (x - fPlus[i] - fPlus[j] + curFValue) /(step * step) ;
        }
    
 
    
    for (int i = 0; i < varCount; i++)
        for (int j = 0; j < i; j++)
            hessian(i,j) = hessian(j, i);
 
 
    return hessian;
}
 
int Hessian()
{
    /*Чтобы вычислить якобиан для другого набора функций, создайте эти функции
     и запишите в массив func указатели на них */
 
    std::ofstream out("output.txt");
    out.setf(std::ios_base::fixed, std::ios_base::floatfield);
    out.precision(15);
    std::cout.setf(std::ios_base::fixed, std::ios_base::floatfield);
    std::cout.precision(6);
    vector<double> point(2, 1.0) ; // Задание точки, в которой вычисляется якобиан
    double step = 1e-4; //Задание шага метода
 
    matrix<double> hessian  = ComputeHessian(&f, point, step);
    matrix<double> ist(2, 2);
    ist(0, 0) = 2;
    ist(1, 0) = 1;
    ist(0, 1) = 1;
    ist(1, 1) = 0.25;
 
    //Вывод результата
    std::cout << "Hessian of given function in point (";
    for (int i = 0; i < point.size() - 1; i++)
        std::cout << point[i] << ", ";
    std::cout << point[ point.size() - 1] << ") equals: " << std::endl; 
    std::cout << hessian;
 
    return 0;
}