diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index ac3b7371a8962969800e44e6849722ae820dc202..faf8db253379a1e291ee33a4ea9e2464229f2f99 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -13,7 +13,7 @@ add_library(4neuro SHARED
         NetConnection/ConnectionWeightIdentity.cpp
         LearningMethods/ParticleSwarm.cpp
         DataSet/DataSet.cpp
-        ErrorFunction/ErrorFunctions.cpp Network/NeuralNetworkSum.cpp Network/NeuralNetworkSum.h)
+        ErrorFunction/ErrorFunctions.cpp Network/NeuralNetworkSum.cpp Network/NeuralNetworkSum.h Solvers/PDESolver.cpp Solvers/PDESolver.h Neuron/NeuronLogistic_d1.cpp Neuron/NeuronLogistic_d1.h Neuron/NeuronLogistic_d2.cpp Neuron/NeuronLogistic_d2.h)
 
 target_link_libraries(4neuro boost_serialization)
 
diff --git a/src/Neuron/Neuron.cpp b/src/Neuron/Neuron.cpp
index 18e9a7ac4137321a675340198572d90ea19b264e..a230bdeb0391cb6f6ac17d4f9b8a12b75382ebdb 100644
--- a/src/Neuron/Neuron.cpp
+++ b/src/Neuron/Neuron.cpp
@@ -124,6 +124,11 @@ void Neuron::set_idx(size_t value) {
     this->neural_net_index = value;
 }
 
+Neuron* IDifferentiable::get_derivative() {
+    return nullptr;
+}
+
+
 //template<class Archive>
 //void Neuron::serialize(Archive & ar, const unsigned int version) {
 //    ar << this->potential;
diff --git a/src/Neuron/Neuron.h b/src/Neuron/Neuron.h
index 78b461f74b60ef0662e44a0998326fea0baed9c3..7379301f97e2aab290bf9c9f7d0b6d4d064c9bf6 100644
--- a/src/Neuron/Neuron.h
+++ b/src/Neuron/Neuron.h
@@ -233,8 +233,9 @@ public:
 
 
 /**
- * Class serving as an interface providing 'activation_function_get_partial_derivative'
- * and 'activation_function_get_derivative' methods.
+ * Class serving as an interface providing 'activation_function_eval_partial_derivative',
+ * 'activation_function_eval_derivative',  'get_partial_derivative' and
+ * 'get_derivative' methods.
  */
 class IDifferentiable {
     /**
@@ -244,13 +245,20 @@ class IDifferentiable {
      * @return Partial derivative of the activation function according to the
      * 'param_idx'-th parameter
      */
-    virtual double activation_function_get_partial_derivative(int param_idx) = 0;
+    virtual double activation_function_eval_partial_derivative(int param_idx) = 0;
 
     /**
      * Calculates the derivative with respect to the argument, ie the 'potential'
      * @return f'(x), where 'f(x)' is the activation function and 'x' = 'potential'
      */
-    virtual double activation_function_get_derivative( ) = 0;
+    virtual double activation_function_eval_derivative( ) = 0;
+
+    /**
+     * Returns a Neuron pointer object with activation function being the partial derivative of
+     * the activation function of this Neuron object with respect to the argument, i.e. 'potential'
+     * @return
+     */
+    virtual Neuron* get_derivative( );
 }; /* end of IDifferentiable class */
 
  #endif /* NEURON_H_ */
\ No newline at end of file
diff --git a/src/Neuron/NeuronLinear.cpp b/src/Neuron/NeuronLinear.cpp
index 2f9acc91a66806ab61ae25b1d219298d76df8d56..e0ae690d9e1fd7c752a243fdf20f3a4c4e2166ba 100644
--- a/src/Neuron/NeuronLinear.cpp
+++ b/src/Neuron/NeuronLinear.cpp
@@ -35,7 +35,7 @@ void NeuronLinear::activate( ) {
 
 }
 
-double NeuronLinear::activation_function_get_partial_derivative(int param_idx) {
+double NeuronLinear::activation_function_eval_partial_derivative(int param_idx) {
 
     if(param_idx == 0){
         return 1.0;
@@ -48,12 +48,21 @@ double NeuronLinear::activation_function_get_partial_derivative(int param_idx) {
     return 0.0;
 }
 
-double NeuronLinear::activation_function_get_derivative( ) {
+double NeuronLinear::activation_function_eval_derivative( ) {
     double b = this->activation_function_parameters[1];
 
     return b;
 }
 
+Neuron* NeuronLinear::get_derivative() {
+    NeuronLinear* output = nullptr;
+
+    //derivative according to 'x'
+    output = new NeuronLinear(this->activation_function_parameters[1], 0.0);
+
+    return output;
+}
+
 //template<class Archive>
 //void NeuronLinear::serialize(Archive & ar, const unsigned int version) {
 //    ar & boost::serialization::base_object<Neuron>(*this);
diff --git a/src/Neuron/NeuronLinear.h b/src/Neuron/NeuronLinear.h
index c74e150600f36d727b0a960c38c114f0aa9886bf..ea3532431b5012d92259ced4aa08285834265a5a 100644
--- a/src/Neuron/NeuronLinear.h
+++ b/src/Neuron/NeuronLinear.h
@@ -45,19 +45,26 @@ public:
 
     /**
      * Calculates the partial derivative of the activation function
-     * f(x) = b*x + a
+     * f(x) = b*x + a at point x
      * @param[in] param_idx Index of the parameter to calculate derivative of
      * @return Partial derivative of the activation function according to the
      * 'param_idx'-th parameter. For 'param_idx'=0 returns 1, for 'param_idx=1'
      * return 'x' and otherwise returns 0.
      */
-    double activation_function_get_partial_derivative(int param_idx) override;
+    double activation_function_eval_partial_derivative(int param_idx) override;
 
     /**
-     * Calculates d/dx of (b*x + a)
+     * Calculates d/dx of (b*x + a) at point x
      * @return b
      */
-    double activation_function_get_derivative( ) override;
+    double activation_function_eval_derivative( ) override;
+
+    /**
+     * Returns a pointer to a Neuron with derivative as its activation function
+     * @return
+     */
+    Neuron* get_derivative() override;
+
 };
 
 
diff --git a/src/Neuron/NeuronLogistic.cpp b/src/Neuron/NeuronLogistic.cpp
index dabcb04722378f9c281f5348b68f34333158edae..21d5bee4e9635a04b83e1a7c83d805570f698f00 100644
--- a/src/Neuron/NeuronLogistic.cpp
+++ b/src/Neuron/NeuronLogistic.cpp
@@ -33,13 +33,14 @@ void NeuronLogistic::activate( ) {
 
 }
 
-double NeuronLogistic::activation_function_get_partial_derivative(int param_idx ) {
+double NeuronLogistic::activation_function_eval_partial_derivative(int param_idx ) {
 
     double a = this->activation_function_parameters[0];
     double b = this->activation_function_parameters[1];
     double x = this->potential;
 
     if(param_idx == 0){
+        //partial derivative according to parameter 'a'
         double ex = std::pow(E, b - x);
 
         double exa= -std::pow(ex + 1.0, -a);
@@ -47,6 +48,7 @@ double NeuronLogistic::activation_function_get_partial_derivative(int param_idx
         return exa * std::log(ex + 1.0);
     }
     else if(param_idx == 1){
+        //partial derivative according to parameter 'b'
         /**
          * TODO
          * Could be write as activation_function_get_derivative() * -1
@@ -61,7 +63,8 @@ double NeuronLogistic::activation_function_get_partial_derivative(int param_idx
 
 }
 
-double NeuronLogistic::activation_function_get_derivative( ) {
+
+double NeuronLogistic::activation_function_eval_derivative( ) {
 
     double a = this->activation_function_parameters[0];
     double b = this->activation_function_parameters[1];
@@ -72,4 +75,15 @@ double NeuronLogistic::activation_function_get_derivative( ) {
 
     return a * ex * ex2;
 
+}
+
+Neuron* NeuronLogistic::get_derivative() {
+    NeuronLogistic_d1 *output = nullptr;
+    double a = this->activation_function_parameters[0];
+    double b = this->activation_function_parameters[1];
+
+    output = new NeuronLogistic_d1(a, b);
+    output->set_potential( this->potential );
+
+    return output;
 }
\ No newline at end of file
diff --git a/src/Neuron/NeuronLogistic.h b/src/Neuron/NeuronLogistic.h
index b96f2f7a39aa3907e0e46b05b937a0209ec5ec51..f8ddd10759ed190221b65a664b125881c3ad164d 100644
--- a/src/Neuron/NeuronLogistic.h
+++ b/src/Neuron/NeuronLogistic.h
@@ -12,6 +12,7 @@
 
 #include <cmath>
 #include "Neuron.h"
+#include "NeuronLogistic_d1.h"
 #include "../constants.h"
 
 class NeuronLogistic:public Neuron, public IDifferentiable {
@@ -36,25 +37,29 @@ public:
     explicit NeuronLogistic(double a = 0.0, double b = 0.0);
 
     /**
-     * Evaluates '(1 + e^(-x + a))^(-b)' and stores the result into the 'state' property
+     * Evaluates '(1 + e^(-x + b))^(-a)' and stores the result into the 'state' property
      */
     void activate( ) override;
 
     /**
      * Calculates the partial derivative of the activation function
-     * f(x) = (1 + e^(-x + a))^(-b)
+     * f(x) = (1 + e^(-x + b))^(-a)
      * @param[in] param_idx Index of the parameter to calculate derivative of
      * @return Partial derivative of the activation function according to the
      * 'param_idx'-th parameter.
      */
-    double activation_function_get_partial_derivative(int param_idx) override;
-
+    double activation_function_eval_partial_derivative(int param_idx) override;
     /**
-     * Calculates d/dx of (1 + e^(-x + a))^(-b)
-     * @return b * e^(a - x) * [e^(a - x) + 1]^(-b)
+     * Calculates d/dx of (1 + e^(-x + b))^(-a)
+     * @return a * e^(b - x) * [e^(b - x) + 1]^(-a)
      */
-    double activation_function_get_derivative( ) override;
+    double activation_function_eval_derivative( ) override;
 
+    /**
+     * Returns a pointer to a Neuron with derivative as its activation function
+     * @return
+     */
+    Neuron* get_derivative() override;
 };
 
 
diff --git a/src/Neuron/NeuronLogistic_d1.cpp b/src/Neuron/NeuronLogistic_d1.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fd5acd3dc3121e19cbdb87b2f2c120a3c47df7dc
--- /dev/null
+++ b/src/Neuron/NeuronLogistic_d1.cpp
@@ -0,0 +1,91 @@
+/**
+ * DESCRIPTION OF THE FILE
+ *
+ * @author Michal KravÄŤenko
+ * @date 22.7.18 -
+ */
+
+#include "NeuronLogistic_d1.h"
+
+Neuron* NeuronLogistic_d1::get_copy( ){
+    NeuronLogistic_d1* output = new NeuronLogistic_d1( this->activation_function_parameters[0],  this->activation_function_parameters[1]);
+
+    return output;
+}
+
+NeuronLogistic_d1::NeuronLogistic_d1(double a, double b) {
+
+    this->n_activation_function_parameters = 2;
+    this->activation_function_parameters = new double[2];
+
+    this->activation_function_parameters[0] = a;
+    this->activation_function_parameters[1] = b;
+
+    this->edges_in = new std::vector<Connection*>(0);
+    this->edges_out = new std::vector<Connection*>(0);
+}
+
+
+void NeuronLogistic_d1::activate( ) {
+    //a*e^(b - x)*(1 + e^(b - x))^(-1 - a)
+
+    double a = this->activation_function_parameters[0];
+    double b = this->activation_function_parameters[1];
+    double x = this->potential;
+
+    double ex = std::pow(E, b - x);
+    this->state = a * ex * std::pow(1.0 + ex, -a - 1.0);
+
+}
+
+double NeuronLogistic_d1::activation_function_eval_partial_derivative(int param_idx ) {
+
+    double a = this->activation_function_parameters[0];
+    double b = this->activation_function_parameters[1];
+    double x = this->potential;
+
+    if(param_idx == 0){
+        //partial derivative according to parameter 'a'
+        //e^(b - x) (1 + e^(b - x))^(-1 - a)[1 - a log(1 + e^(b - x))]
+        double ex = std::pow(E, b - x);
+        double ex2= std::pow(1.0 + ex, -1.0 - a);
+
+        return ex * ex2 * (1.0 - a * std::log(1.0 + ex));
+    }
+    else if(param_idx == 1){
+        //partial derivative according to parameter 'b'
+        //[(-1 - a) e^(b-x) (1 + e^(b - x))^(-1) + 1] * a e^(b - x) (1 + e^(b - x))^(-1 - a)
+        double ex = std::pow(E, b - x);
+        double ex2 = std::pow(ex + 1.0, -a - 1.0);
+        double ex3 = 1.0 / (1.0 + ex);
+
+        return a * ex * ex2 * (1.0 - (a + 1.0) * ex * ex3);
+    }
+
+    return 0.0;
+}
+
+double NeuronLogistic_d1::activation_function_eval_derivative() {
+    //[(1 + a) e^(b-x) (1 + e^(b - x))^(-1) - 1]a e^(b - x) (1 + e^(b - x))^(-1 - a)
+
+    double a = this->activation_function_parameters[0];
+    double b = this->activation_function_parameters[1];
+    double x = this->potential;
+
+    double ex = std::pow(E, b - x);
+    double ex2 = std::pow(ex + 1.0, -a - 1.0);
+    double ex3 = 1.0 / (1.0 + ex);
+
+    return -a * ex * ex2 * (1.0 - (a + 1.0) * ex * ex3);
+}
+
+ Neuron* NeuronLogistic_d1::get_derivative() {
+    //[(1 + a) e^(b-x) (1 + e^(b - x))^(-1) - 1]a e^(b - x) (1 + e^(b - x))^(-1 - a)
+    NeuronLogistic_d2* output = nullptr;
+
+     double a = this->activation_function_parameters[0];
+     double b = this->activation_function_parameters[1];
+     output = new NeuronLogistic_d2(a, b);
+
+    return output;
+}
\ No newline at end of file
diff --git a/src/Neuron/NeuronLogistic_d1.h b/src/Neuron/NeuronLogistic_d1.h
new file mode 100644
index 0000000000000000000000000000000000000000..66e5ed8b9615cc92b3fec08aef8e295e85a02330
--- /dev/null
+++ b/src/Neuron/NeuronLogistic_d1.h
@@ -0,0 +1,64 @@
+/**
+ * DESCRIPTION OF THE FILE
+ *
+ * @author Michal KravÄŤenko
+ * @date 22.7.18 -
+ */
+
+#ifndef INC_4NEURO_NEURONLOGISTIC_D1_H
+#define INC_4NEURO_NEURONLOGISTIC_D1_H
+
+#include <cmath>
+#include "Neuron.h"
+#include "NeuronLogistic_d2.h"
+#include "../constants.h"
+
+class NeuronLogistic_d1:public Neuron, public IDifferentiable {
+    friend class boost::serialization::access;
+protected:
+    template<class Archive>
+    void serialize(Archive & ar, const unsigned int version){
+        //TODO separate implementation to NeuronLogistic_d1.cpp!
+        ar & boost::serialization::base_object<Neuron>(*this);
+    };
+public:
+
+    Neuron* get_copy( );
+
+    /**
+     * Constructs the object of the Logistic neuron with activation function
+     * f(x) = a*e^(b - x)*(1 + e^(b - x))^(-1 - a)
+     * @param[in] a First coefficient, stored in activation_function_parameters[0]
+     * @param[in] b Second coefficient, stored in activation_function_parameters[1]
+     */
+    explicit NeuronLogistic_d1(double a = 0.0, double b = 0.0);
+
+    /**
+     * Evaluates 'a*e^(b - x)*(1 + e^(b - x))^(-1 - a)' and stores the result into the 'state' property
+     */
+    void activate( ) override;
+
+    /**
+     * Calculates the partial derivative of the activation function
+     * f(x) = a*e^(b - x)*(1 + e^(b - x))^(-1 - a)
+     * @param[in] param_idx Index of the parameter to calculate derivative of
+     * @return Partial derivative of the activation function according to the
+     * 'param_idx'-th parameter.
+     */
+    double activation_function_eval_partial_derivative(int param_idx) override;
+
+    /**
+     * Calculates d/dx of  [a*e^(b - x)*(1 + e^(b - x))^(-1 - a)]
+     * @return [[(1 + a) e^(b-x) (1 + e^(b - x))^(-1) - 1]a e^(b - x) (1 + e^(b - x))^(-1 - a)]
+     */
+    double activation_function_eval_derivative( ) override;
+
+    /**
+     * Returns a pointer to a Neuron with derivative as its activation function
+     * @return
+     */
+    Neuron* get_derivative() override;
+};
+
+
+#endif //INC_4NEURO_NEURONLOGISTIC_D1_H
diff --git a/src/Neuron/NeuronLogistic_d2.cpp b/src/Neuron/NeuronLogistic_d2.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bfa30f086566cbcc7d2da234c7956abd8c8d687d
--- /dev/null
+++ b/src/Neuron/NeuronLogistic_d2.cpp
@@ -0,0 +1,87 @@
+/**
+ * DESCRIPTION OF THE FILE
+ *
+ * @author Michal KravÄŤenko
+ * @date 22.7.18 -
+ */
+
+#include "NeuronLogistic_d2.h"
+
+Neuron* NeuronLogistic_d2::get_copy( ){
+    NeuronLogistic_d2* output = new NeuronLogistic_d2( this->activation_function_parameters[0],  this->activation_function_parameters[1]);
+
+    return output;
+}
+
+NeuronLogistic_d2::NeuronLogistic_d2(double a, double b) {
+
+    this->n_activation_function_parameters = 2;
+    this->activation_function_parameters = new double[2];
+
+    this->activation_function_parameters[0] = a;
+    this->activation_function_parameters[1] = b;
+
+    this->edges_in = new std::vector<Connection*>(0);
+    this->edges_out = new std::vector<Connection*>(0);
+}
+
+void NeuronLogistic_d2::activate( ) {
+    //[(1 + a) e^(b-x) (1 + e^(b - x))^(-1) - 1]a e^(b - x) (1 + e^(b - x))^(-1 - a)
+
+    double a = this->activation_function_parameters[0];
+    double b = this->activation_function_parameters[1];
+    double x = this->potential;
+
+    double ex = std::pow(E, b - x);
+    double ex2 = std::pow(ex + 1.0, -a - 1.0);
+    double ex3 = 1.0 / (1.0 + ex);
+
+    this->state = -a * ex * ex2 * (1.0 - (a + 1.0) * ex * ex3);
+
+}
+
+double NeuronLogistic_d2::activation_function_eval_partial_derivative(int param_idx ) {
+
+    double a = this->activation_function_parameters[0];
+    double b = this->activation_function_parameters[1];
+    double x = this->potential;
+
+    if(param_idx == 0){
+        //partial derivative according to parameter 'a'
+        //(e^b (e^(b - x) + 1)^(-a) (a^2 (-e^b) log(e^(b - x) + 1) + a e^x log(e^(b - x) + 1) + 2 a e^b - e^x))/(e^b + e^x)^2
+        double eb = std::pow(E, b);
+        double ex = std::pow(E, x);
+        double ebx= eb / ex;
+        double ebxa = std::pow(ebx + 1.0, -a);
+        double lbx = std::log(ebx + 1.0);
+
+        return eb * ebxa * (a * lbx * ( a * (-eb) + ex ) + 2.0 * a * eb - ex) / ((eb + ex) * (eb + ex));
+    }
+    else if(param_idx == 1){
+        //partial derivative according to parameter 'b'
+        //-(a e^b (e^(b - x) + 1)^(-a) (a^2 e^(2 b) - 3 a e^(b + x) - e^(b + x) + e^(2 x)))/(e^b + e^x)^3
+
+        double eb = std::pow(E, b);
+        double ex = std::pow(E, x);
+        double ebx= eb / ex;
+        double ebxa = std::pow(ebx + 1.0, -a);
+
+        return  - a * eb * ebxa * (a * a * eb * eb - 3.0 * a * ebx - ebx + ex * ex) / ((eb + ex) * (eb + ex) * (eb + ex));
+    }
+
+    return 0.0;
+}
+
+double NeuronLogistic_d2::activation_function_eval_derivative() {
+    //(a e^b (1 + e^(b - x))^(-a) (a^2 e^(2 b) + e^(2 x) - e^(b + x) - 3 a e^(b + x)))/(e^b + e^x)^3
+    double a = this->activation_function_parameters[0];
+    double b = this->activation_function_parameters[1];
+    double x = this->potential;
+
+    double eb = std::pow(E, b);
+    double ex = std::pow(E, x);
+    double ebx= eb / ex;
+    double ebxa = std::pow(ebx + 1.0, -a);
+
+    return  a * eb * ebxa * (a * a * eb * eb - 3.0 * a * ebx - ebx + ex * ex) / ((eb + ex) * (eb + ex) * (eb + ex));
+}
diff --git a/src/Neuron/NeuronLogistic_d2.h b/src/Neuron/NeuronLogistic_d2.h
new file mode 100644
index 0000000000000000000000000000000000000000..d40bce0bda0425473f113b9fc4e2b1f419160be6
--- /dev/null
+++ b/src/Neuron/NeuronLogistic_d2.h
@@ -0,0 +1,59 @@
+/**
+ * DESCRIPTION OF THE FILE
+ *
+ * @author Michal KravÄŤenko
+ * @date 22.7.18 -
+ */
+
+#ifndef INC_4NEURO_NEURONLOGISTIC_D2_H
+#define INC_4NEURO_NEURONLOGISTIC_D2_H
+
+
+#include <cmath>
+#include "Neuron.h"
+#include "../constants.h"
+
+class NeuronLogistic_d2:public Neuron, public IDifferentiable {
+    friend class boost::serialization::access;
+protected:
+    template<class Archive>
+    void serialize(Archive & ar, const unsigned int version){
+        //TODO separate implementation to NeuronLogistic_d1.cpp!
+        ar & boost::serialization::base_object<Neuron>(*this);
+    };
+public:
+
+    Neuron* get_copy( );
+
+    /**
+     * Constructs the object of the Logistic neuron with activation function
+     * f(x) = [(1 + a) e^(b-x) (1 + e^(b - x))^(-1) - 1]a e^(b - x) (1 + e^(b - x))^(-1 - a)
+     * @param[in] a First coefficient, stored in activation_function_parameters[0]
+     * @param[in] b Second coefficient, stored in activation_function_parameters[1]
+     */
+    explicit NeuronLogistic_d2(double a = 0.0, double b = 0.0);
+
+    /**
+     * Evaluates 'a*e^(b - x)*(1 + e^(b - x))^(-1 - a)' and stores the result into the 'state' property
+     */
+    void activate( ) override;
+
+    /**
+     * Calculates the partial derivative of the activation function
+     * f(x) = [(1 + a) e^(b-x) (1 + e^(b - x))^(-1) - 1]a e^(b - x) (1 + e^(b - x))^(-1 - a)
+     * @param[in] param_idx Index of the parameter to calculate derivative of
+     * @return Partial derivative of the activation function according to the
+     * 'param_idx'-th parameter.
+     */
+    double activation_function_eval_partial_derivative(int param_idx) override;
+
+    /**
+     * Calculates d/dx of  [(1 + a) e^(b-x) (1 + e^(b - x))^(-1) - 1]a e^(b - x) (1 + e^(b - x))^(-1 - a)
+     * @return (a e^b (1 + e^(b - x))^(-a) (a^2 e^(2 b) + e^(2 x) - e^(b + x) - 3 a e^(b + x)))/(e^b + e^x)^3
+     */
+    double activation_function_eval_derivative( ) override;
+
+};
+
+
+#endif //INC_4NEURO_NEURONLOGISTIC_D2_H
diff --git a/src/Neuron/NeuronTanh.cpp b/src/Neuron/NeuronTanh.cpp
index 80071c33a60925cceddc6e7859d6d7a9020204fd..cef2a07faedf044e19f06c368b0729d5c1756bb0 100644
--- a/src/Neuron/NeuronTanh.cpp
+++ b/src/Neuron/NeuronTanh.cpp
@@ -32,7 +32,7 @@ void NeuronTanh::activate( ) {
 
 }
 
-double NeuronTanh::activation_function_get_partial_derivative(int param_idx) {
+double NeuronTanh::activation_function_eval_partial_derivative(int param_idx) {
 
     double a = this->activation_function_parameters[0];
     double x = this->potential;
@@ -53,7 +53,7 @@ double NeuronTanh::activation_function_get_partial_derivative(int param_idx) {
 
 }
 
-double NeuronTanh::activation_function_get_derivative( ) {
+double NeuronTanh::activation_function_eval_derivative( ) {
 
     double a = this->activation_function_parameters[0];
     double x = this->potential;
@@ -61,4 +61,8 @@ double NeuronTanh::activation_function_get_derivative( ) {
     double exi = std::pow(E, 2.0 * a) + std::pow(E, 2.0 * x);
 
     return ex / (exi * exi);
+}
+
+Neuron* NeuronTanh::get_derivative() {
+    return nullptr;
 }
\ No newline at end of file
diff --git a/src/Neuron/NeuronTanh.h b/src/Neuron/NeuronTanh.h
index a1e9b91fc58534be09a14ce61d537dfe77e8c4ce..3d5ce6d8aabf4ad6bcabb55193612606dde7241d 100644
--- a/src/Neuron/NeuronTanh.h
+++ b/src/Neuron/NeuronTanh.h
@@ -47,7 +47,7 @@ public:
      * @return Partial derivative of the activation function according to the
      * 'param_idx'-th parameter.
      */
-    double activation_function_get_partial_derivative(int param_idx) override;
+    double activation_function_eval_partial_derivative(int param_idx) override;
 
     /**
      * TODO
@@ -55,7 +55,13 @@ public:
      * Calculates d/dx of (e^(x-a) - e^(a-x))/(e^(x-a) + e^(a-x))
      * @return a * e^(b - x) * [e^(b - x) + 1]^(-a)
      */
-    double activation_function_get_derivative( ) override;
+    double activation_function_eval_derivative( ) override;
+
+    /**
+     * Returns a pointer to a Neuron with derivative as its activation function
+     * @return
+     */
+    Neuron* get_derivative() override;
 };
 
 
diff --git a/src/Solvers/PDESolver.cpp b/src/Solvers/PDESolver.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..86030da681bff15f7f979de020f5b03d94d075b1
--- /dev/null
+++ b/src/Solvers/PDESolver.cpp
@@ -0,0 +1,8 @@
+/**
+ * DESCRIPTION OF THE FILE
+ *
+ * @author Michal KravÄŤenko
+ * @date 22.7.18 - 
+ */
+
+#include "PDESolver.h"
diff --git a/src/Solvers/PDESolver.h b/src/Solvers/PDESolver.h
new file mode 100644
index 0000000000000000000000000000000000000000..74c2968036d45048191b8c83f19fb72ddf527ddc
--- /dev/null
+++ b/src/Solvers/PDESolver.h
@@ -0,0 +1,81 @@
+/**
+ * File containing methods for quick & simple formulation of PDEs as a system of Neural Networks with shared weights
+ *
+ * @author Michal KravÄŤenko
+ * @date 22.7.18 -
+ */
+
+#ifndef INC_4NEURO_PDESOLVER_H
+#define INC_4NEURO_PDESOLVER_H
+
+
+#include "../DataSet/DataSet.h"
+#include "../Network/NeuralNetwork.h"
+
+/**
+ * class representing a multi-index of partial derivatives
+ */
+struct MultiIndex{
+    /**
+     * total number of variables
+     */
+    unsigned int dim;
+
+    /**
+     * a vector containing degrees of partial derivatives with respect to each variable
+     */
+    std::vector<unsigned int> partial_derivatives_degrees;
+};
+
+class PDESolver {
+
+public:
+    /**
+     * default constructor
+     */
+    PDESolver();
+
+    /**
+     * default destructor
+     */
+    ~PDESolver();
+
+    /**
+     * adds a new boundary condition in the form of
+     * [d^{k}/d^{alpha}{x_{alpha_1}x_{alpha_2}...x_{alpha_k}} y(x_1, x_2, ..., x_n) = f(x_1, x_2, ..., x_n)]
+     * to the system.
+     * @param alpha
+     * @param conditions
+     */
+    void add_boundary_condition(MultiIndex &alpha, DataSet &conditions);
+
+
+    /**
+     * adds a new partial derivative in the form of
+     * [d^{k}/d^{alpha}{x_{alpha_1}x_{alpha_2}...x_{alpha_k}} y(x_1, x_2, ..., x_n)] * beta
+     * to the PDE.
+     * @param alpha
+     * @param beta
+     */
+    void add_to_equation(MultiIndex &alpha, double beta);
+
+    /**
+     * sets the right hand side so PDE(x_1, x_2, ..., x_n) = rhs(x_1, x_2, ..., x_n)
+     * @param rhs
+     */
+    void set_quation_rhs(DataSet &rhs);
+
+    /**
+     * solves the PDE with its boundary conditions via the particle swarm algorithm
+     */
+    void solve_particle_swarm( );
+
+    /**
+     * returns the pointer to the object representing the solution y(x_1, x_2, ..., x_n)
+     * @return
+     */
+    NeuralNetwork* get_solution( );
+};
+
+
+#endif //INC_4NEURO_PDESOLVER_H