diff --git a/src/DataSet/DataSet.cpp b/src/DataSet/DataSet.cpp
index 661ec79ed0d543bdd0d0e3577afb6fcc3995dbd3..dac9f2f0f03e8c0e1be92e154af7efd739fdb525 100644
--- a/src/DataSet/DataSet.cpp
+++ b/src/DataSet/DataSet.cpp
@@ -40,7 +40,7 @@ DataSet::DataSet(double lower_bound, double upper_bound, unsigned int size, doub
     this->add_isotropic_data(lower_bound, upper_bound, size, output);
 }
 
-DataSet::DataSet(std::vector<double> bounds, unsigned int no_elems_in_one_dim, std::vector<double> (*output_func)(std::vector<double>), unsigned int output_dim) {
+DataSet::DataSet(std::vector<double> &bounds, unsigned int no_elems_in_one_dim, std::vector<double> (*output_func)(std::vector<double>&), unsigned int output_dim) {
     std::vector<std::pair<std::vector<double>, std::vector<double>>> new_data_vec;
     this->data = new_data_vec;
     this->input_dim = bounds.size()/2;
@@ -51,7 +51,7 @@ DataSet::DataSet(std::vector<double> bounds, unsigned int no_elems_in_one_dim, s
 }
 
 
-void DataSet::add_data_pair(std::vector<double> inputs, std::vector<double> outputs) {
+void DataSet::add_data_pair(std::vector<double> &inputs, std::vector<double> &outputs) {
     if(inputs.size() != this->input_dim) {
         throw InvalidDimension("Bad input dimension.");
     } else if(outputs.size() != this->output_dim) {
@@ -81,7 +81,7 @@ void DataSet::add_isotropic_data(double lower_bound, double upper_bound, unsigne
     this->n_elements += size;
 }
 
-void DataSet::add_isotropic_data(std::vector<double> bounds, unsigned int no_elems_in_one_dim, std::vector<double> (*output_func)(std::vector<double>)) {
+void DataSet::add_isotropic_data(std::vector<double> &bounds, unsigned int no_elems_in_one_dim, std::vector<double> (*output_func)(std::vector<double>&)) {
     // TODO add check of dataset dimensions
 
     std::vector<std::vector<double>> grid;
@@ -100,7 +100,7 @@ void DataSet::add_isotropic_data(std::vector<double> bounds, unsigned int no_ele
 
     grid = this->cartesian_product(&grid);
 
-    for(const auto& vec : grid) {
+    for(auto vec : grid) {
         this->n_elements++;
         this->data.emplace_back(std::make_pair(vec, output_func(vec)));
     }
@@ -114,11 +114,11 @@ size_t DataSet::get_n_elements() {
     return this->n_elements;
 }
 
-unsigned int DataSet::get_input_dim() {
+size_t DataSet::get_input_dim() {
     return this->input_dim;
 }
 
-unsigned int DataSet::get_output_dim() {
+size_t DataSet::get_output_dim() {
     return this->output_dim;
 }
 
@@ -142,7 +142,7 @@ void DataSet::print_data() {
     }
 }
 
-void DataSet::store_text(std::string file_path) {
+void DataSet::store_text(std::string &file_path) {
     //TODO check if stream was successfully opened
     std::ofstream ofs(file_path);
     boost::archive::text_oarchive oa(ofs);
diff --git a/src/DataSet/DataSet.h b/src/DataSet/DataSet.h
index 2d47988d4415b657d3c3497b8d393939071b2ea1..11654263cfd3704a4486cd054487c91c8a033a89 100644
--- a/src/DataSet/DataSet.h
+++ b/src/DataSet/DataSet.h
@@ -49,12 +49,12 @@ private:
     /**
      * Dimension of the input
      */
-    unsigned int input_dim = 0;
+    size_t input_dim = 0;
 
     /**
      * Dimension of the output
      */
-    unsigned int output_dim = 0;
+    size_t output_dim = 0;
 
     /**
      * Stored data in the format of pairs of corresponding
@@ -169,7 +169,14 @@ public:
      */
     DataSet(double lower_bound, double upper_bound, unsigned int size, double output);
 
-    DataSet(std::vector<double> bounds, unsigned int no_elems_in_one_dim, std::vector<double> (*output_func)(std::vector<double>), unsigned int output_dim);
+    /**
+     *
+     * @param bounds
+     * @param no_elems_in_one_dim
+     * @param output_func
+     * @param output_dim
+     */
+    DataSet(std::vector<double> &bounds, unsigned int no_elems_in_one_dim, std::vector<double> (*output_func)(std::vector<double>&), unsigned int output_dim);
 
     /**
      * Getter for number of elements
@@ -181,14 +188,14 @@ public:
      * Returns the input dimension
      * @return Input dimension
      */
-    unsigned int get_input_dim();
+    size_t get_input_dim();
 
 
     /**
      * Return the output dimension
      * @return Output dimension
      */
-    unsigned int get_output_dim();
+    size_t get_output_dim();
 
     /**
      * Getter for the data structure
@@ -201,7 +208,7 @@ public:
      * @param inputs Vector of input data
      * @param outputs Vector of output data corresponding to the input data
      */
-    void add_data_pair(std::vector<double> inputs, std::vector<double> outputs);
+    void add_data_pair(std::vector<double> &inputs, std::vector<double> &outputs);
 
     //TODO expand method to generate multiple data types - chebyshev etc.
     /**
@@ -230,7 +237,7 @@ public:
      * @param size Number of input-output pairs generated
      * @param output_func Function determining output value
      */
-    void add_isotropic_data(std::vector<double> bounds, unsigned int no_elems_in_one_dim, std::vector<double> (*output_func)(std::vector<double>));
+    void add_isotropic_data(std::vector<double> &bounds, unsigned int no_elems_in_one_dim, std::vector<double> (*output_func)(std::vector<double>&));
 
     //TODO Chebyshev - ch. interpolation points, i-th point = cos(i*alpha) from 0 to pi
 
@@ -242,7 +249,7 @@ public:
     /**
      * Stores the DataSet object to the binary file
      */
-    void store_text(std::string file_path);
+    void store_text(std::string &file_path);
 };
 
 #endif //INC_4NEURO_DATASET_H
diff --git a/src/LearningMethods/ParticleSwarm.cpp b/src/LearningMethods/ParticleSwarm.cpp
index 3f1caefc95336e6e71a8f283a94dcafdc7022ceb..c21365c8749ba0f696ac640c3aacd031e662b817 100644
--- a/src/LearningMethods/ParticleSwarm.cpp
+++ b/src/LearningMethods/ParticleSwarm.cpp
@@ -131,8 +131,8 @@ double Particle::change_coordinate(double w, double c1, double c2, std::vector<d
     for(unsigned int i = 0; i < this->coordinate_dim; ++i) {
         vel_mem = w * (*this->velocity)[i]
                   + c1 * this->r1 * ((*this->optimal_coordinate)[i] - (*this->coordinate)[i])
-                  + c2 * this->r2 * (glob_min_coord[i] - (*this->coordinate)[i]);
-//                  + (c1+c2)/2 * this->r3 * ((*random_global_best)[i] - (*this->coordinate)[i]);
+                  + c2 * this->r2 * (glob_min_coord[i] - (*this->coordinate)[i])
+                  + (c1+c2)/2 * this->r3 * ((*random_global_best)[i] - (*this->coordinate)[i]);
 
         if ((*this->coordinate)[i] + vel_mem > this->domain_bounds[2 * i + 1]) {
             this->randomize_velocity();
diff --git a/src/Solvers/DESolver.cpp b/src/Solvers/DESolver.cpp
index dfc9b2aa38e1b86092bd96b932db89ff7a2c1e9c..061bc2a21e19e5668d45c3b245702fa539c8169a 100644
--- a/src/Solvers/DESolver.cpp
+++ b/src/Solvers/DESolver.cpp
@@ -361,10 +361,10 @@ NeuralNetwork* DESolver::get_solution() {
     return this->solution;
 }
 
-double DESolver::eval_equation( size_t equation_idx, std::vector<double> &weight_and_biases, std::vector<double> &input ) {
+double DESolver::eval_equation( size_t equation_idx, std::vector<double> *weight_and_biases, std::vector<double> &input ) {
     std::vector<double> output(1);
 
-    this->differential_equations->at( equation_idx )->eval_single( input, output, &weight_and_biases );
+    this->differential_equations->at( equation_idx )->eval_single( input, output, weight_and_biases );
 
 //    printf("Input: ");
 //    for( auto e: input ){
diff --git a/src/Solvers/DESolver.h b/src/Solvers/DESolver.h
index e9e953605f642bff8b89c4411bad5734ed45df05..0ae33b3afda2cc0adc0ab5a9b039909007f1b44d 100644
--- a/src/Solvers/DESolver.h
+++ b/src/Solvers/DESolver.h
@@ -165,7 +165,7 @@ public:
     /**
      * For testing purposes only
      */
-     double eval_equation( size_t equation_idx, std::vector<double> &weights_and_biases, std::vector<double> &input );
+     double eval_equation( size_t equation_idx, std::vector<double> *weights_and_biases, std::vector<double> &input );
 
      /**
       * For testing purposes only
diff --git a/src/examples/net_test_ode_1.cpp b/src/examples/net_test_ode_1.cpp
index 106ee9d1a9912c40d1f69ec0624197d6ab811c79..b140db5e238b34eead1076bc6ca8cfc0580ba420 100644
--- a/src/examples/net_test_ode_1.cpp
+++ b/src/examples/net_test_ode_1.cpp
@@ -541,7 +541,7 @@ void test_analytical_gradient_y(std::vector<double> &guess, double accuracy, siz
     delete conjugate_direction_prev;
 }
 
-void test_odr(double accuracy, size_t n_inner_neurons, size_t train_size, double ds, double de, size_t n_test_points, double ts, double te, size_t max_iters, size_t n_particles){
+void test_ode(double accuracy, size_t n_inner_neurons, size_t train_size, double ds, double de, size_t n_test_points, double ts, double te, size_t max_iters, size_t n_particles){
 
     /* SOLVER SETUP */
     size_t n_inputs = 1;
@@ -664,8 +664,6 @@ void test_odr(double accuracy, size_t n_inner_neurons, size_t train_size, double
         printf("\tRepresentation test %6d, error of eq1: %10.8f, error of eq2: %10.8f, error of eq3: %10.8f, total error: %10.8f\n", (int)testi, std::sqrt(test_error_eq1), std::sqrt(test_error_eq2), std::sqrt(test_error_eq3), (total_error_analytical - total_error_de_solver) * (total_error_analytical - total_error_de_solver));
     }
 
-return;
-
     /* PARTICLE SWARM TRAINING METHOD SETUP */
 
     //must encapsulate each of the partial error functions
@@ -675,11 +673,11 @@ return;
         domain_bounds[2 * i + 1] = 10.0;
     }
 
-    double c1 = 0.5, c2 = 0.8, w = 0.7;
+    double c1 = 1.7, c2 = 1.7, w = 0.7;
 
 
 
-    double gamma = 0.5, epsilon = 0.02, delta = 0.9;
+    double gamma = 0.5, epsilon = 0.02, delta = 1.0;
 
     solver_01.solve_via_particle_swarm( domain_bounds, c1, c2, w, n_particles, max_iters, gamma, epsilon, delta );
 
@@ -712,7 +710,7 @@ return;
 int main() {
 
     unsigned int n_inner_neurons = 2;
-    unsigned int train_size = 150;
+    unsigned int train_size = 10;
     double accuracy = 1e-4;
     double ds = 0.0;
     double de = 4.0;
@@ -722,8 +720,8 @@ int main() {
     double te = de + 2;
 
     size_t particle_swarm_max_iters = 1000;
-    size_t n_particles = 10;
-    test_odr(accuracy, n_inner_neurons, train_size, ds, de, test_size, ts, te, particle_swarm_max_iters, n_particles);
+    size_t n_particles = 100;
+    test_ode(accuracy, n_inner_neurons, train_size, ds, de, test_size, ts, te, particle_swarm_max_iters, n_particles);
 
 //    bool optimize_weights = true;
 //    bool optimize_biases = true;
diff --git a/src/examples/net_test_pde_1.cpp b/src/examples/net_test_pde_1.cpp
index 998cca0430f61085c0853d47ac5b9b435e58c92b..d2c1ffebfa65e8558164d9391e99032afa298829 100644
--- a/src/examples/net_test_pde_1.cpp
+++ b/src/examples/net_test_pde_1.cpp
@@ -17,21 +17,15 @@
 
 #include <random>
 #include <iostream>
+#include <fstream>
 
 #include "../../include/4neuro.h"
 #include "../Solvers/DESolver.h"
 
 
-void test_odr(size_t n_inner_neurons){
+void test_pde(double accuracy, size_t n_inner_neurons, size_t train_size, double ds, double de, size_t n_test_points, double ts, double te, size_t max_iters, size_t n_particles){
 
     /* solution properties */
-    size_t train_size = 10;
-    double d1_s = 0.0, d1_e = 1.0;
-
-    /* swarm optimizer properties */
-    unsigned int max_iters = 100;
-    unsigned int n_particles = 10;
-
 
     /* do not change below */
     size_t n_inputs = 2;
@@ -42,9 +36,13 @@ void test_odr(size_t n_inner_neurons){
     MultiIndex alpha_00( n_inputs );
     MultiIndex alpha_01( n_inputs );
     MultiIndex alpha_20( n_inputs );
+
     alpha_00.set_partial_derivative(0, 0);
-    alpha_01.set_partial_derivative(0, 1);
-    alpha_20.set_partial_derivative(2, 0);
+    alpha_00.set_partial_derivative(1, 0);
+
+    alpha_01.set_partial_derivative(1, 1);
+
+    alpha_20.set_partial_derivative(0, 2);
 
     /* the governing differential equation */
     solver_01.add_to_differential_equation( 0, alpha_20,  1.0 );
@@ -61,44 +59,36 @@ void test_odr(size_t n_inner_neurons){
     double frac, x, t;
 
     /* TRAIN DATA FOR THE GOVERNING DE */
-    std::vector<std::pair<std::vector<double>, std::vector<double>>> data_vec_g;
+    std::vector<double> test_bounds_2d = {ds, de, ds, de};
 
+    /* GOVERNING EQUATION RHS */
+    auto f1 = [](std::vector<double>&input) -> std::vector<double> {
+        std::vector<double> output(1);
+        output[0] = 0.0;
+        return output;
+    };
+    DataSet ds_00(test_bounds_2d, train_size, f1, 1);
 
-    /* ISOTROPIC TRAIN SET */
-    frac = (d1_e - d1_s) / (train_size - 1);
-    for(size_t i = 0; i < train_size; ++i){
-        inp = {0.0, 0.0};
-        out = {0.0};
-        data_vec_g.emplace_back(std::make_pair(inp, out));
-    }
-    DataSet ds_00(&data_vec_g);
-
-    /* ISOTROPIC TRAIN DATA FOR THE FIRST DIRICHLET BC */
     std::vector<std::pair<std::vector<double>, std::vector<double>>> data_vec_t;
-    frac = (d1_e - d1_s) / (train_size - 1);
-    for(size_t i = 0; i < train_size; ++i){
-        t = i * frac;
-        inp = { 0.0, t };
-        out = { std::sin( t ) };
-
-        data_vec_t.emplace_back(std::make_pair(inp, out));
-    }
-    DataSet ds_t(&data_vec_t);
-
-    /* ISOTROPIC TRAIN DATA FOR THE SECOND DIRICHLET BC */
     std::vector<std::pair<std::vector<double>, std::vector<double>>> data_vec_x;
-    frac = (d1_e - d1_s) / (train_size - 1);
-    for(size_t i = 0; i < train_size; ++i){
-        x = i * frac;
-        inp = { x, 0.0 };
-        out = { std::pow(E, -0.707106781 * x) * std::sin( -0.707106781 * x )};
+    /* ISOTROPIC TRAIN SET */
+    frac = (de - ds) / (train_size - 1);
+    for(unsigned int i = 0; i < train_size; ++i){
+        inp = {0.0, frac * i};
+        out = {std::sin(inp[1])};
+        data_vec_t.emplace_back(std::make_pair(inp, out));
 
+        inp = {frac * i, 0.0};
+        out = {std::pow(E, -0.707106781 * inp[0]) * std::sin( -0.707106781 * inp[0] )};
         data_vec_x.emplace_back(std::make_pair(inp, out));
+
     }
+    DataSet ds_t(&data_vec_t);
     DataSet ds_x(&data_vec_x);
 
 
 
+
     /* Placing the conditions into the solver */
     solver_01.set_error_function( 0, ErrorFunctionType::ErrorFuncMSE, &ds_00 );
     solver_01.set_error_function( 1, ErrorFunctionType::ErrorFuncMSE, &ds_t );
@@ -112,40 +102,74 @@ void test_odr(size_t n_inner_neurons){
     //must encapsulate each of the partial error functions
     double *domain_bounds = new double[ 6 * n_inner_neurons ];
     for(unsigned int i = 0; i < 3 * n_inner_neurons; ++i){
-        domain_bounds[2 * i] = -800.0;
-        domain_bounds[2 * i + 1] = 800.0;
+        domain_bounds[2 * i] = -20.0;
+        domain_bounds[2 * i + 1] = 20.0;
     }
 
-    double c1 = 0.5, c2 = 0.75, w = 0.8;
+    double c1 = 1.7, c2 = 1.7, w = 0.7;
+    double gamma = 0.5, epsilon = 0.02, delta = 0.9;
 
+    solver_01.solve_via_particle_swarm( domain_bounds, c1, c2, w, n_particles, max_iters, gamma, epsilon, delta );
+    /* PRACTICAL END OF THE EXAMPLE */
 
 
-    double gamma = 0.5, epsilon = 0.02, delta = 0.9;
 
-    solver_01.solve_via_particle_swarm( domain_bounds, c1, c2, w, n_particles, max_iters, gamma, epsilon, delta );
+    /* SOLUTION EXPORT */
+    printf("Exporting solution & error files...");
 
-//    NeuralNetwork *solution = solver_01.get_solution();
-//    std::vector<double> parameters(3 * n_inner_neurons);//w1, a1, b1, w2, a2, b2, ... , wm, am, bm
-//    std::vector<double> *weight_params = solution->get_parameter_ptr_weights();
-//    std::vector<double> *biases_params = solution->get_parameter_ptr_biases();
-//    for(size_t i = 0; i < n_inner_neurons; ++i){
-//        parameters[3 * i] = weight_params->at(i);
-//        parameters[3 * i + 1] = weight_params->at(i + n_inner_neurons);
-//        parameters[3 * i + 2] = biases_params->at(i);
-//    }
-//
-//    unsigned int n_test_points = 150;
-//    std::vector<double> input(1), output(1);
-//    double x;
-//    for(unsigned int i = 0; i < n_test_points; ++i){
-//
-//        x = i * ((d1_e - d1_s) / (n_test_points - 1)) + d1_s;
-//        input[0] = x;
-//
-//        solution->eval_single(input, output);
-//
-//        std::cout << i + 1 << " " << x << " " << std::pow(E, -2*x) * (3*x + 1)<< " " << output[0] << " " << std::pow(E, -2*x) * (1 - 6*x)<< " " << eval_approx_df(x, n_inner_neurons, parameters) << " " << 4 * std::pow(E, -2*x) * (3*x - 2)<< " " << eval_approx_ddf(x, n_inner_neurons, parameters) << std::endl;
-//    }
+    NeuralNetwork *solution = solver_01.get_solution();
+    std::vector<double> *weight_params = solution->get_parameter_ptr_weights();
+    std::vector<double> *biases_params = solution->get_parameter_ptr_biases();
+
+    /* solution itself */
+    DataSet test_set_1(test_bounds_2d, n_test_points, f1, 1);
+
+    std::vector<double> input, output(1);
+    std::ofstream ofs("data_2d_pde1_y.txt", std::ofstream::out);
+    for(auto tp: *test_set_1.get_data()){
+        input = tp.first;
+
+        solution->eval_single(input, output);
+
+        ofs << input[0] << " " << input[1] << " " << output[0] << std::endl;
+    }
+    ofs.close();
+
+    /* governing equation error */
+    ofs = std::ofstream("data_2d_pde1_first_equation_error.txt", std::ofstream::out);
+    for(auto tp: *test_set_1.get_data()){
+        input = tp.first;
+
+        double eq_value = solver_01.eval_equation(0, nullptr, input);
+
+        ofs << input[0] << " " << input[1] << " " << std::fabs(eq_value) << std::endl;
+    }
+    ofs.close();
+
+    /* ISOTROPIC TEST SET FOR BOUNDARY CONDITIONS */
+    frac = (de - ds) / (n_test_points - 1);
+    /* first boundary condition & its error */
+    ofs = std::ofstream("data_1d_pde1_yt.txt", std::ofstream::out);
+    std::ofstream ofs2("data_1d_pde1_yx.txt", std::ofstream::out);
+    for(unsigned int i = 0; i < n_test_points; ++i){
+        double x = frac * i;
+        double t = frac * i;
+
+        double yt = std::sin(t);
+        double yx = std::pow(E, -0.707106781 * x) * std::sin( -0.707106781 * x );
+
+        input = {0.0, t};
+        solution->eval_single( input, output, nullptr );
+        ofs << i + 1 << " " << t << " " << yt << " " << output[0] << " " << std::fabs(output[0] - yt) << std::endl;
+
+        input = {x, 0.0};
+        solution->eval_single( input, output, nullptr );
+        ofs2 << i + 1 << " " << x << " " << yx << " " << output[0] << " " << std::fabs(output[0] - yx) << std::endl;
+    }
+    ofs2.close();
+    ofs.close();
+
+    printf("done!\n");
 
 
     delete [] domain_bounds;
@@ -153,9 +177,19 @@ void test_odr(size_t n_inner_neurons){
 
 int main() {
 
-    unsigned int n_inner_neurons = 2;
+    unsigned int n_inner_neurons = 6;
+    unsigned int train_size = 20;
+    double accuracy = 1e-4;
+    double ds = 0.0;
+    double de = 1.0;
+
+    unsigned int test_size = 100;
+    double ts = ds;
+    double te = de;
 
-    test_odr(n_inner_neurons);
+    size_t particle_swarm_max_iters = 1000;
+    size_t n_particles = 100;
+    test_pde(accuracy, n_inner_neurons, train_size, ds, de, test_size, ts, te, particle_swarm_max_iters, n_particles);
 
 
     return 0;