#include <iostream>
#include <octave/oct.h>
#include <octave/builtin-defun-decls.h>

int main (void) {
    
    Matrix A(4,4);
    
    for (octave_idx_type i = 0; i < 4; i++)
       for (octave_idx_type j = 0; j < 4; j++)
        A(i,j) = 1.0 / (static_cast<double> (i) +static_cast<double> ( j ) + 1.0 ) ;
    
        ColumnVector b(4,1.0);
    
    // solve Ax=b
    
    octave_value_list in;
    in_mat(0)=A;
    
    octave_value_list in1;
    in_vec(0) =b;
    
    //octave_value_list norm =Fnorm(in);
    
    octave_value_list sol=Fmldivide(in_mat,in_vec); // ???? don't know what are the right auguments here
    
   // double norm_of_the_matrix = sol(0).double_value ();
  //  double norm_of_the_matrix = norm(0).double_value ();
      //  ColumnVector x = Fmldivide(octave_value(b),1)(0).column_vector_value();
    //ColumnVector x = do_binary_op(octave_value::op_ldiv,A,b,1).column_vector_value();
        std::cout << "A = " << std::endl << A<< std::endl
                  << "b = " << std::endl << b << std::endl
                  << "sol = " << std::endl << sol << std::endl;
                  //<< "norm_of_the_matrix = " << std::endl <<norm_of_the_matrix << std::endl;
        return 0;
}