diff --git a/Makefile b/Makefile
index 0735553..633e5bc 100644
--- a/Makefile
+++ b/Makefile
@@ -7,7 +7,7 @@ CFLAGS += -O3 -flto=auto -march=native
 CXXFLAGS = $(CFLAGS)
 LDFLAGS = $(CXXFLAGS)
 LDLIBS  = -lm
-OBJS = partdiff.o askparams.o
+OBJS = partdiff.o argument_parser.o calculation_arguments.o tensor.o
 
 default: all
 
diff --git a/argument_parser.cpp b/argument_parser.cpp
new file mode 100644
index 0000000..b5907e5
--- /dev/null
+++ b/argument_parser.cpp
@@ -0,0 +1,194 @@
+#include "argument_parser.hpp"
+#include "calculation_options.hpp"
+#include <format>
+#include <print>
+#include <sstream>
+
+namespace partdiff {
+
+  template <typename T, typename U = std::underlying_type_t<T>>
+  U to_underlying(T v) {
+    return static_cast<U>(v);
+  }
+
+  static constexpr uint64_t max_interlines = 10240;
+  static constexpr uint64_t max_iteration = 200000;
+  static constexpr uint64_t max_threads = 1024;
+  static constexpr double min_accuracy = 1e-4;
+  static constexpr double max_accuracy = 1e-20;
+
+  using calculation_method = calculation_options::calculation_method;
+  using perturbation_function = calculation_options::perturbation_function;
+  using termination_condition = calculation_options::termination_condition;
+
+  argument_parser::argument_parser(const int argc, char const *argv[])
+    : app_name(argv[0]),
+      args(argv + 1, argv + argc) {
+    this->fill_argument_descriptions();
+    this->ask_params();
+  }
+
+  calculation_options argument_parser::get_options() {
+    return this->options;
+  }
+
+  argument_parser::argument_description argument_parser::get_description(std::size_t index) const {
+    return this->argument_descriptions[index];
+  }
+
+  argument_parser::argument_description argument_parser::get_description(argument_index index) const {
+    return this->get_description(to_underlying(index));
+  }
+
+  void argument_parser::usage() const {
+    const auto get_name = [](const std::string &input) { return std::format("  - {:11}", input + ":"); };
+
+    std::print("Usage: {}", this->app_name);
+    for (std::size_t i = 0; i <= to_underlying(argument_index::term_dummy); i++) {
+      std::print(" [{}]", get_description(i).name);
+    }
+    std::println("");
+    std::println("");
+    for (std::size_t i = 0; i <= to_underlying(argument_index::term_dummy); i++) {
+      const std::string description = this->get_description(i).description_for_usage.value_or("<invalid>");
+      std::println("{}{}", get_name(this->get_description(i).name), description);
+    }
+    std::println("Example: {} 1 2 100 1 2 100", app_name);
+  }
+
+  void argument_parser::ask_params() {
+    if (this->args.size() < 6) {
+      usage();
+      exit(EXIT_SUCCESS);
+    }
+    for (std::size_t i = 0; i <= to_underlying(argument_index::termination); i++) {
+      parse_param(i, args[i]);
+    }
+    if (this->options.termination == termination_condition::accuracy) {
+      parse_param(argument_index::term_accuracy, args[5]);
+      this->options.term_iteration = partdiff::max_iteration;
+    } else {
+      parse_param(argument_index::term_iteration, args[5]);
+      this->options.term_accuracy = 0.0;
+    }
+  }
+
+  void argument_parser::parse_param(argument_index index, std::string &input) {
+    this->parse_param(to_underlying(index), input);
+  }
+
+  void argument_parser::parse_param(std::size_t index, std::string &input) {
+    if (!this->get_description(index).read_from_string(input)) {
+      this->usage();
+      exit(EXIT_FAILURE);
+    }
+  }
+
+  void argument_parser::fill_argument_descriptions() {
+
+    auto scientific_double = [](double val) {
+      auto temp = std::format("{:.0e}", val);
+      int epos = temp.find("e");
+      std::string mantissa_str = temp.substr(0, epos);
+      std::string exponent_str = temp.substr(epos + 1, temp.length() - epos - 1);
+      int exponent = stoi(exponent_str);
+      return mantissa_str + "e" + std::to_string(exponent);
+    };
+
+    constexpr int indent_width = 17;
+    const std::string indent = std::format("{:{}s}", "", indent_width);
+
+    auto number = &(this->options.number);
+    this->add_argument_description("num", number, std::format("number of threads (1 .. {:d})", partdiff::max_threads),
+                                   [number] { return (*number >= 1 && *number <= partdiff::max_threads); });
+
+    auto method = &(this->options.method);
+    this->add_argument_description(
+        "method", method,
+        std::format("calculation method (1 .. 2)\n"
+                    "{0}{1:d}: Gauß-Seidel\n"
+                    "{0}{2:d}: Jacobi",
+                    indent, to_underlying(calculation_method::gauss_seidel), to_underlying(calculation_method::jacobi)),
+        [method] { return (*method == calculation_method::gauss_seidel || *method == calculation_method::jacobi); });
+
+    auto interlines = &(this->options.interlines);
+    this->add_argument_description("lines", interlines,
+                                   std::format("number of interlines (0 .. {1:d})\n"
+                                               "{0}matrixsize = (interlines * 8) + 9",
+                                               indent, partdiff::max_interlines),
+                                   [interlines] { return (*interlines <= partdiff::max_interlines); });
+
+    auto pert_func = &(this->options.pert_func);
+    this->add_argument_description(
+        "func", pert_func,
+        std::format("perturbation function (1 .. 2)\n"
+                    "{0}{1:d}: f(x,y) = 0\n"
+                    "{0}{2:d}: f(x,y) = 2 * pi^2 * sin(pi * x) * sin(pi * y)",
+                    indent, to_underlying(perturbation_function::f0), to_underlying(perturbation_function::fpisin)),
+        [pert_func] {
+          return (*pert_func == perturbation_function::f0 || *pert_func == perturbation_function::fpisin);
+        });
+
+    auto termination = &(this->options.termination);
+    this->add_argument_description("term", termination,
+                                   std::format("termination condition ( 1.. 2)\n"
+                                               "{0}{1:d}: sufficient accuracy\n"
+                                               "{0}{2:d}: number of iterations",
+                                               indent, to_underlying(termination_condition::accuracy),
+                                               to_underlying(termination_condition::iterations)),
+                                   [termination] {
+                                     return (*termination == termination_condition::accuracy ||
+                                             *termination == termination_condition::iterations);
+                                   });
+
+    this->add_argument_description("acc/iter",
+                                   std::format("depending on term:\n"
+                                               "{0}accuracy:  {1:s} .. {2:s}\n"
+                                               "{0}iterations:    1 .. {3:d}\n",
+                                               indent, scientific_double(partdiff::min_accuracy),
+                                               scientific_double(partdiff::max_accuracy), partdiff::max_iteration));
+
+    auto term_accuracy = &(this->options.term_accuracy);
+    this->add_argument_description("acc", term_accuracy, std::nullopt, [term_accuracy] {
+      return (*term_accuracy >= partdiff::max_accuracy && *term_accuracy <= partdiff::min_accuracy);
+    });
+
+    auto term_iteration = &(this->options.term_iteration);
+    this->add_argument_description("iter", term_iteration, std::nullopt, [term_iteration] {
+      return (*term_iteration >= 1 && *term_iteration <= partdiff::max_iteration);
+    });
+  }
+
+  void argument_parser::add_argument_description(std::string name, std::optional<std::string> description_for_usage) {
+    argument_description arg_desc;
+    arg_desc.name = name;
+    arg_desc.description_for_usage = description_for_usage;
+    this->argument_descriptions.push_back(arg_desc);
+  }
+
+  template <class T>
+  void argument_parser::add_argument_description(std::string name, T *target,
+                                                 std::optional<std::string> description_for_usage,
+                                                 std::function<bool()> check) {
+    argument_description arg_desc;
+    arg_desc.name = name;
+    arg_desc.target = target;
+    arg_desc.read_from_string = [target = arg_desc.target, check](const std::string &input) {
+      T *casted_ptr = std::any_cast<T *>(target);
+      bool valid_input = false;
+      std::istringstream iss(input);
+      if constexpr (std::is_enum_v<T>) {
+        std::underlying_type_t<T> temp;
+        valid_input = static_cast<bool>(iss >> temp);
+        *casted_ptr = static_cast<T>(temp);
+      } else {
+        valid_input = static_cast<bool>(iss >> *casted_ptr);
+      }
+      valid_input &= check();
+      return valid_input;
+    };
+    arg_desc.description_for_usage = description_for_usage;
+    this->argument_descriptions.push_back(arg_desc);
+  }
+
+} // namespace partdiff
diff --git a/argument_parser.hpp b/argument_parser.hpp
new file mode 100644
index 0000000..ac6a7b9
--- /dev/null
+++ b/argument_parser.hpp
@@ -0,0 +1,53 @@
+#include "calculation_options.hpp"
+#include <any>
+#include <functional>
+#include <optional>
+#include <string>
+#include <vector>
+
+namespace partdiff {
+
+  class argument_parser {
+    public:
+    argument_parser(const int argc, char const *argv[]);
+    calculation_options get_options();
+
+    private:
+    struct argument_description {
+      std::any target;
+      std::string name;
+      std::optional<std::string> description_for_usage;
+      std::function<bool(const std::string &input)> read_from_string = [](auto) { return false; };
+    };
+
+    enum class argument_index : std::size_t {
+      number = 0,
+      method = 1,
+      interlines = 2,
+      pert_func = 3,
+      termination = 4,
+      term_dummy = 5,
+      term_accuracy = 6,
+      term_iteration = 7
+    };
+
+    calculation_options options;
+    std::string app_name;
+    std::vector<std::string> args;
+    std::vector<argument_description> argument_descriptions;
+    argument_description get_description(std::size_t index) const;
+    argument_description get_description(argument_index index) const;
+    void usage() const;
+    void ask_params();
+    void parse_param(std::size_t index, std::string &input);
+    void parse_param(argument_index index, std::string &input);
+    void ask_param(std::size_t index);
+    void ask_param(argument_index index);
+    void fill_argument_descriptions();
+    template <class T>
+    void add_argument_description(std::string name, T *target, std::optional<std::string> description_for_usage,
+                                  std::function<bool()> check);
+    void add_argument_description(std::string name, std::optional<std::string> description_for_usage);
+  };
+
+} // namespace partdiff
diff --git a/askparams.cpp b/askparams.cpp
deleted file mode 100644
index f2b3e9c..0000000
--- a/askparams.cpp
+++ /dev/null
@@ -1,184 +0,0 @@
-#include "partdiff.h"
-
-namespace partdiff {
-
-  namespace askparams {
-
-    using calculation_method = calculation_options::calculation_method;
-    using perturbation_function = calculation_options::perturbation_function;
-    using termination_condition = calculation_options::termination_condition;
-
-    argument_parser::argument_parser(const int argc, char const *argv[])
-      : app_name(argv[0]),
-        args(argv + 1, argv + argc) {
-      this->fill_argument_descriptions();
-      this->ask_params();
-    }
-
-    calculation_options argument_parser::get_options() {
-      return this->options;
-    }
-
-    argument_parser::argument_description argument_parser::get_description(std::size_t index) const {
-      return this->argument_descriptions[index];
-    }
-
-    argument_parser::argument_description argument_parser::get_description(argument_index index) const {
-      return this->get_description(to_underlying(index));
-    }
-
-    void argument_parser::usage() const {
-      const auto get_name = [](const std::string &input) { return std::format("  - {:11}", input + ":"); };
-
-      std::print("Usage: {}", this->app_name);
-      for (std::size_t i = 0; i <= to_underlying(argument_index::term_dummy); i++) {
-        std::print(" [{}]", get_description(i).name);
-      }
-      std::println("");
-      std::println("");
-      for (std::size_t i = 0; i <= to_underlying(argument_index::term_dummy); i++) {
-        const std::string description = this->get_description(i).description_for_usage.value_or("<invalid>");
-        std::println("{}{}", get_name(this->get_description(i).name), description);
-      }
-      std::println("Example: {} 1 2 100 1 2 100", app_name);
-    }
-
-    void argument_parser::ask_params() {
-      if (this->args.size() < 6) {
-        usage();
-        exit(EXIT_SUCCESS);
-      }
-      for (std::size_t i = 0; i <= to_underlying(argument_index::termination); i++) {
-        parse_param(i, args[i]);
-      }
-      if (this->options.termination == termination_condition::accuracy) {
-        parse_param(argument_index::term_accuracy, args[5]);
-        this->options.term_iteration = partdiff::max_iteration;
-      } else {
-        parse_param(argument_index::term_iteration, args[5]);
-        this->options.term_accuracy = 0.0;
-      }
-    }
-
-    void argument_parser::parse_param(argument_index index, std::string &input) {
-      this->parse_param(to_underlying(index), input);
-    }
-
-    void argument_parser::parse_param(std::size_t index, std::string &input) {
-      if (!this->get_description(index).read_from_string(input)) {
-        this->usage();
-        exit(EXIT_FAILURE);
-      }
-    }
-
-    void argument_parser::fill_argument_descriptions() {
-
-      auto scientific_double = [](double val) {
-        auto temp = std::format("{:.0e}", val);
-        int epos = temp.find("e");
-        std::string mantissa_str = temp.substr(0, epos);
-        std::string exponent_str = temp.substr(epos + 1, temp.length() - epos - 1);
-        int exponent = stoi(exponent_str);
-        return mantissa_str + "e" + std::to_string(exponent);
-      };
-
-      constexpr int indent_width = 17;
-      const std::string indent = std::format("{:{}s}", "", indent_width);
-
-      auto number = &(this->options.number);
-      this->add_argument_description("num", number, std::format("number of threads (1 .. {:d})", partdiff::max_threads),
-                                     [number] { return (*number >= 1 && *number <= partdiff::max_threads); });
-
-      auto method = &(this->options.method);
-      this->add_argument_description(
-          "method", method,
-          std::format("calculation method (1 .. 2)\n"
-                      "{0}{1:d}: Gauß-Seidel\n"
-                      "{0}{2:d}: Jacobi",
-                      indent, to_underlying(calculation_method::gauss_seidel),
-                      to_underlying(calculation_method::jacobi)),
-          [method] { return (*method == calculation_method::gauss_seidel || *method == calculation_method::jacobi); });
-
-      auto interlines = &(this->options.interlines);
-      this->add_argument_description("lines", interlines,
-                                     std::format("number of interlines (0 .. {1:d})\n"
-                                                 "{0}matrixsize = (interlines * 8) + 9",
-                                                 indent, partdiff::max_interlines),
-                                     [interlines] { return (*interlines <= partdiff::max_interlines); });
-
-      auto pert_func = &(this->options.pert_func);
-      this->add_argument_description(
-          "func", pert_func,
-          std::format("perturbation function (1 .. 2)\n"
-                      "{0}{1:d}: f(x,y) = 0\n"
-                      "{0}{2:d}: f(x,y) = 2 * pi^2 * sin(pi * x) * sin(pi * y)",
-                      indent, to_underlying(perturbation_function::f0), to_underlying(perturbation_function::fpisin)),
-          [pert_func] {
-            return (*pert_func == perturbation_function::f0 || *pert_func == perturbation_function::fpisin);
-          });
-
-      auto termination = &(this->options.termination);
-      this->add_argument_description("term", termination,
-                                     std::format("termination condition ( 1.. 2)\n"
-                                                 "{0}{1:d}: sufficient accuracy\n"
-                                                 "{0}{2:d}: number of iterations",
-                                                 indent, to_underlying(termination_condition::accuracy),
-                                                 to_underlying(termination_condition::iterations)),
-                                     [termination] {
-                                       return (*termination == termination_condition::accuracy ||
-                                               *termination == termination_condition::iterations);
-                                     });
-
-      this->add_argument_description("acc/iter",
-                                     std::format("depending on term:\n"
-                                                 "{0}accuracy:  {1:s} .. {2:s}\n"
-                                                 "{0}iterations:    1 .. {3:d}\n",
-                                                 indent, scientific_double(partdiff::min_accuracy),
-                                                 scientific_double(partdiff::max_accuracy), partdiff::max_iteration));
-
-      auto term_accuracy = &(this->options.term_accuracy);
-      this->add_argument_description("acc", term_accuracy, std::nullopt, [term_accuracy] {
-        return (*term_accuracy >= partdiff::max_accuracy && *term_accuracy <= partdiff::min_accuracy);
-      });
-
-      auto term_iteration = &(this->options.term_iteration);
-      this->add_argument_description("iter", term_iteration, std::nullopt, [term_iteration] {
-        return (*term_iteration >= 1 && *term_iteration <= partdiff::max_iteration);
-      });
-    }
-
-    void argument_parser::add_argument_description(std::string name, std::optional<std::string> description_for_usage) {
-      argument_description arg_desc;
-      arg_desc.name = name;
-      arg_desc.description_for_usage = description_for_usage;
-      this->argument_descriptions.push_back(arg_desc);
-    }
-
-    template <class T>
-    void argument_parser::add_argument_description(std::string name, T *target,
-                                                   std::optional<std::string> description_for_usage,
-                                                   std::function<bool()> check) {
-      argument_description arg_desc;
-      arg_desc.name = name;
-      arg_desc.target = target;
-      arg_desc.read_from_string = [target = arg_desc.target, check](const std::string &input) {
-        T *casted_ptr = std::any_cast<T *>(target);
-        bool valid_input = false;
-        std::istringstream iss(input);
-        if constexpr (std::is_enum_v<T>) {
-          std::underlying_type_t<T> temp;
-          valid_input = static_cast<bool>(iss >> temp);
-          *casted_ptr = static_cast<T>(temp);
-        } else {
-          valid_input = static_cast<bool>(iss >> *casted_ptr);
-        }
-        valid_input &= check();
-        return valid_input;
-      };
-      arg_desc.description_for_usage = description_for_usage;
-      this->argument_descriptions.push_back(arg_desc);
-    }
-
-  } // namespace askparams
-
-} // namespace partdiff
diff --git a/calculation_arguments.cpp b/calculation_arguments.cpp
new file mode 100644
index 0000000..80237f1
--- /dev/null
+++ b/calculation_arguments.cpp
@@ -0,0 +1,41 @@
+#include "calculation_arguments.hpp"
+#include "calculation_options.hpp"
+
+namespace partdiff {
+
+  using calculation_options = calculation_options;
+  using calculation_method = calculation_options::calculation_method;
+  using perturbation_function = calculation_options::perturbation_function;
+  using termination_condition = calculation_options::termination_condition;
+
+  calculation_arguments::calculation_arguments(const calculation_options &options) : pert_func(options.pert_func) {
+    this->N = (options.interlines * 8) + 9 - 1;
+    this->num_matrices = (options.method == calculation_method::jacobi) ? 2 : 1;
+    this->h = 1.0 / this->N;
+    this->matrices = tensor(num_matrices, N + 1, N + 1);
+    this->init_matrices();
+  }
+
+  void calculation_arguments::init_matrices() {
+    for (uint64_t g = 0; g < this->num_matrices; g++) {
+      for (uint64_t i = 0; i <= N; i++) {
+        for (uint64_t j = 0; j <= N; j++) {
+          this->matrices[g, i, j] = 0.0;
+        }
+      }
+    }
+    if (this->pert_func == perturbation_function::f0) {
+      for (uint64_t g = 0; g < this->num_matrices; g++) {
+        for (uint64_t i = 0; i <= N; i++) {
+          this->matrices[g, i, 0] = 1.0 - (h * i);
+          this->matrices[g, i, N] = h * i;
+          this->matrices[g, 0, i] = 1.0 - (h * i);
+          this->matrices[g, N, i] = h * i;
+        }
+        this->matrices[g, N, 0] = 0.0;
+        this->matrices[g, 0, N] = 0.0;
+      }
+    }
+  }
+
+} // namespace partdiff
diff --git a/calculation_arguments.hpp b/calculation_arguments.hpp
new file mode 100644
index 0000000..e071f72
--- /dev/null
+++ b/calculation_arguments.hpp
@@ -0,0 +1,21 @@
+#pragma once
+
+#include "calculation_options.hpp"
+#include "tensor.hpp"
+
+namespace partdiff {
+
+  struct calculation_arguments {
+
+    uint64_t N;
+    uint64_t num_matrices;
+    double h;
+    tensor matrices;
+    calculation_arguments(const calculation_options &);
+
+    private:
+    calculation_options::perturbation_function pert_func;
+    void init_matrices();
+  };
+
+} // namespace partdiff
diff --git a/calculation_options.hpp b/calculation_options.hpp
new file mode 100644
index 0000000..df8c945
--- /dev/null
+++ b/calculation_options.hpp
@@ -0,0 +1,20 @@
+#pragma once
+
+#include <cstdint>
+
+namespace partdiff {
+
+  struct calculation_options {
+    enum class calculation_method : uint64_t { gauss_seidel = 1, jacobi = 2 };
+    enum class perturbation_function : uint64_t { f0 = 1, fpisin = 2 };
+    enum class termination_condition : uint64_t { accuracy = 1, iterations = 2 };
+    uint64_t number;
+    uint64_t interlines;
+    calculation_method method;
+    perturbation_function pert_func;
+    termination_condition termination;
+    uint64_t term_iteration;
+    double term_accuracy;
+  };
+
+} // namespace partdiff
diff --git a/calculation_results.hpp b/calculation_results.hpp
new file mode 100644
index 0000000..f4a3299
--- /dev/null
+++ b/calculation_results.hpp
@@ -0,0 +1,16 @@
+#pragma once
+
+#include <chrono>
+
+namespace partdiff {
+
+  struct calculation_results {
+    using time_point = std::chrono::time_point<std::chrono::high_resolution_clock>;
+    int m;
+    uint64_t stat_iteration;
+    double stat_accuracy;
+    time_point start_time;
+    time_point end_time;
+  };
+
+} // namespace partdiff
diff --git a/justfile b/justfile
index da39b62..3556657 100644
--- a/justfile
+++ b/justfile
@@ -5,4 +5,4 @@ build:
     make
 
 format:
-    clang-format --style=file -i *.h *.cpp
+    clang-format --style=file -i *.cpp *.hpp
diff --git a/partdiff.cpp b/partdiff.cpp
index f438878..fce03bc 100644
--- a/partdiff.cpp
+++ b/partdiff.cpp
@@ -1,114 +1,30 @@
+#include "argument_parser.hpp"
+#include "calculation_arguments.hpp"
+#include "calculation_options.hpp"
+#include "calculation_results.hpp"
 #include <cmath>
-#include <numbers>
-#include <utility>
-
-#include "partdiff.h"
+#include <format>
+#include <print>
 
 namespace partdiff {
 
-  using argument_parser = askparams::argument_parser;
-  using calculation_options = askparams::calculation_options;
+  using argument_parser = argument_parser;
+  using calculation_options = calculation_options;
   using calculation_method = calculation_options::calculation_method;
   using perturbation_function = calculation_options::perturbation_function;
   using termination_condition = calculation_options::termination_condition;
-  using tensor = calculation_arguments::tensor;
 
   static constexpr double pi = std::numbers::pi;
   static constexpr double two_pi_square = (2 * pi * pi);
 
-  inline tensor::tensor(std::size_t num_matrices, std::size_t num_rows, std::size_t num_cols)
-    : num_matrices(num_matrices),
-      num_rows(num_rows),
-      num_cols(num_cols) {
-    auto size = num_matrices * num_rows * num_cols;
-    try {
-      data = new double[size];
-    } catch (std::bad_alloc &) {
-      auto size_bytes = size * sizeof(double);
-      std::println("Memory failure! (Requested {} bytes)", size_bytes);
-      exit(EXIT_FAILURE);
-    }
-  }
-
-  inline tensor::tensor(const tensor &other)
-    : num_matrices(other.num_matrices),
-      num_rows(other.num_rows),
-      num_cols(other.num_cols),
-      data(other.data) {}
-
-  inline tensor::tensor(tensor &&other) noexcept
-    : num_matrices(other.num_matrices),
-      num_rows(other.num_rows),
-      num_cols(other.num_cols),
-      data(std::exchange(other.data, nullptr)) {}
-
-  inline tensor &tensor::operator=(const tensor &other) {
-    return *this = tensor(other);
-  }
-
-  inline tensor &tensor::operator=(tensor &&other) noexcept // move assignment
-  {
-    std::swap(data, other.data);
-    num_matrices = other.num_matrices;
-    num_cols = other.num_cols;
-    num_rows = other.num_rows;
-    return *this;
-  }
+  static calculation_results calculate(calculation_arguments &arguments, const calculation_options &options) {
 
-  inline tensor::~tensor() {
-    if (data) {
-      delete[] data;
-      data = nullptr;
-    }
-  }
+    const auto now = std::chrono::high_resolution_clock::now;
 
-  inline double &tensor::operator()(std::size_t matrix, std::size_t row, std::size_t col) {
-    return data[(num_cols * num_rows * matrix) + (num_cols * row) + (col)];
-  }
+    const auto start_time = now();
 
-  inline double tensor::operator()(std::size_t matrix, std::size_t row, std::size_t col) const {
-    return data[(num_cols * num_rows * matrix) + (num_cols * row) + (col)];
-  }
-
-  calculation_arguments::calculation_arguments(const calculation_options &options) : pert_func(options.pert_func) {
-    this->N = (options.interlines * 8) + 9 - 1;
-    this->num_matrices = (options.method == calculation_method::jacobi) ? 2 : 1;
-    this->h = 1.0 / this->N;
-    this->matrices = tensor(num_matrices, N + 1, N + 1);
-    this->init_matrices();
-  }
-
-  void calculation_arguments::init_matrices() {
-    for (uint64_t g = 0; g < this->num_matrices; g++) {
-      for (uint64_t i = 0; i <= N; i++) {
-        for (uint64_t j = 0; j <= N; j++) {
-          this->matrices(g, i, j) = 0.0;
-        }
-      }
-    }
-    if (this->pert_func == perturbation_function::f0) {
-      for (uint64_t g = 0; g < this->num_matrices; g++) {
-        for (uint64_t i = 0; i <= N; i++) {
-          this->matrices(g, i, 0) = 1.0 - (h * i);
-          this->matrices(g, i, N) = h * i;
-          this->matrices(g, 0, i) = 1.0 - (h * i);
-          this->matrices(g, N, i) = h * i;
-        }
-        this->matrices(g, N, 0) = 0.0;
-        this->matrices(g, 0, N) = 0.0;
-      }
-    }
-  }
-
-  calculation_results::calculation_results() {
-    this->m = 0;
-    this->stat_iteration = 0;
-    this->stat_accuracy = 0;
-  }
-
-  static void calculate(calculation_arguments &arguments, calculation_results &results,
-                        const calculation_options &options) {
-    results.start_time = std::chrono::high_resolution_clock::now();
+    uint64_t stat_iteration = 0;
+    double stat_accuracy = 0.0;
 
     const int N = arguments.N;
     const double h = arguments.h;
@@ -140,25 +56,25 @@ namespace partdiff {
         }
 
         for (int j = 1; j < N; j++) {
-          double star = 0.25 * (arguments.matrices(m2, i - 1, j) + arguments.matrices(m2, i, j - 1) +
-                                arguments.matrices(m2, i, j + 1) + arguments.matrices(m2, i + 1, j));
+          double star = 0.25 * (arguments.matrices[m2, i - 1, j] + arguments.matrices[m2, i, j - 1] +
+                                arguments.matrices[m2, i, j + 1] + arguments.matrices[m2, i + 1, j]);
 
           if (options.pert_func == perturbation_function::fpisin) {
             star += fpisin_i * std::sin(pih * (double)j);
           }
 
           if (options.termination == termination_condition::accuracy || term_iteration == 1) {
-            double residuum = arguments.matrices(m2, i, j) - star;
+            double residuum = arguments.matrices[m2, i, j] - star;
             residuum = std::fabs(residuum);
             maxresiduum = std::max(residuum, maxresiduum);
           }
 
-          arguments.matrices(m1, i, j) = star;
+          arguments.matrices[m1, i, j] = star;
         }
       }
 
-      results.stat_iteration++;
-      results.stat_accuracy = maxresiduum;
+      stat_iteration++;
+      stat_accuracy = maxresiduum;
 
       const int temp = m1;
       m1 = m2;
@@ -173,8 +89,10 @@ namespace partdiff {
       }
     }
 
-    results.m = m2;
-    results.end_time = std::chrono::high_resolution_clock::now();
+    const auto end_time = now();
+
+    calculation_results results = {m2, stat_iteration, stat_accuracy, start_time, end_time};
+    return results;
   }
 
   static void display_statistics(const calculation_arguments &arguments, const calculation_results &results,
@@ -208,7 +126,7 @@ namespace partdiff {
     for (int y = 0; y < 9; y++) {
       for (int x = 0; x < 9; x++) {
         std::print(" {:.4f}",
-                   arguments.matrices(results.m, y * (options.interlines + 1), x * (options.interlines + 1)));
+                   arguments.matrices[results.m, y * (options.interlines + 1), x * (options.interlines + 1)]);
       }
       std::println("");
     }
@@ -218,16 +136,15 @@ namespace partdiff {
 
 using calculation_arguments = partdiff::calculation_arguments;
 using calculation_results = partdiff::calculation_results;
-using argument_parser = partdiff::askparams::argument_parser;
-using calculation_options = partdiff::askparams::calculation_options;
+using argument_parser = partdiff::argument_parser;
+using calculation_options = partdiff::calculation_options;
 
 int main(const int argc, char const *argv[]) {
   argument_parser parser(argc, argv);
   calculation_options options = parser.get_options();
   calculation_arguments arguments(options);
-  calculation_results results;
 
-  calculate(arguments, results, options);
+  calculation_results results = calculate(arguments, options);
 
   display_statistics(arguments, results, options);
   display_matrix(arguments, results, options);
diff --git a/partdiff.h b/partdiff.h
deleted file mode 100644
index e5f7154..0000000
--- a/partdiff.h
+++ /dev/null
@@ -1,127 +0,0 @@
-#pragma once
-
-#include <any>
-#include <chrono>
-#include <format>
-#include <functional>
-#include <iomanip>
-#include <iostream>
-#include <optional>
-#include <print>
-#include <sstream>
-#include <vector>
-
-namespace partdiff {
-
-  static constexpr uint64_t max_interlines = 10240;
-  static constexpr uint64_t max_iteration = 200000;
-  static constexpr uint64_t max_threads = 1024;
-  static constexpr double min_accuracy = 1e-4;
-  static constexpr double max_accuracy = 1e-20;
-
-  template <typename T, typename U = std::underlying_type_t<T>>
-  U to_underlying(T v) {
-    return static_cast<U>(v);
-  }
-
-  namespace askparams {
-
-    struct calculation_options {
-      enum class calculation_method : uint64_t { gauss_seidel = 1, jacobi = 2 };
-      enum class perturbation_function : uint64_t { f0 = 1, fpisin = 2 };
-      enum class termination_condition : uint64_t { accuracy = 1, iterations = 2 };
-      uint64_t number;
-      uint64_t interlines;
-      calculation_method method;
-      perturbation_function pert_func;
-      termination_condition termination;
-      uint64_t term_iteration;
-      double term_accuracy;
-    };
-
-    class argument_parser {
-      public:
-      argument_parser(const int argc, char const *argv[]);
-      calculation_options get_options();
-
-      private:
-      struct argument_description {
-        std::any target;
-        std::string name;
-        std::optional<std::string> description_for_usage;
-        std::function<bool(const std::string &input)> read_from_string = [](auto) { return false; };
-      };
-
-      enum class argument_index : std::size_t {
-        number = 0,
-        method = 1,
-        interlines = 2,
-        pert_func = 3,
-        termination = 4,
-        term_dummy = 5,
-        term_accuracy = 6,
-        term_iteration = 7
-      };
-
-      calculation_options options;
-      std::string app_name;
-      std::vector<std::string> args;
-      std::vector<argument_description> argument_descriptions;
-      argument_description get_description(std::size_t index) const;
-      argument_description get_description(argument_index index) const;
-      void usage() const;
-      void ask_params();
-      void parse_param(std::size_t index, std::string &input);
-      void parse_param(argument_index index, std::string &input);
-      void ask_param(std::size_t index);
-      void ask_param(argument_index index);
-      void fill_argument_descriptions();
-      template <class T>
-      void add_argument_description(std::string name, T *target, std::optional<std::string> description_for_usage,
-                                    std::function<bool()> check);
-      void add_argument_description(std::string name, std::optional<std::string> description_for_usage);
-    };
-
-  } // namespace askparams
-
-  struct calculation_arguments {
-
-    class tensor {
-      public:
-      tensor() {};
-      tensor(std::size_t num_matrices, std::size_t num_rows, std::size_t num_cols);
-      tensor(const tensor &other);
-      tensor(tensor &&other) noexcept;
-      tensor &operator=(const tensor &other);
-      tensor &operator=(tensor &&other) noexcept;
-      ~tensor();
-      double &operator()(std::size_t matrix, std::size_t row, std::size_t col);
-      double operator()(std::size_t matrix, std::size_t row, std::size_t col) const;
-
-      private:
-      std::size_t num_matrices, num_rows, num_cols;
-      double *data = nullptr;
-    };
-
-    uint64_t N;
-    uint64_t num_matrices;
-    double h;
-    tensor matrices;
-    calculation_arguments(const askparams::calculation_options &);
-
-    private:
-    askparams::calculation_options::perturbation_function pert_func;
-    void init_matrices();
-  };
-
-  struct calculation_results {
-    using time_point = std::chrono::time_point<std::chrono::high_resolution_clock>;
-    uint64_t m;
-    uint64_t stat_iteration;
-    double stat_accuracy;
-    time_point start_time;
-    time_point end_time;
-    calculation_results();
-  };
-
-} // namespace partdiff
diff --git a/tensor.cpp b/tensor.cpp
new file mode 100644
index 0000000..4db8308
--- /dev/null
+++ b/tensor.cpp
@@ -0,0 +1,61 @@
+#include "tensor.hpp"
+#include <print>
+#include <utility>
+
+namespace partdiff {
+
+  tensor::tensor(std::size_t num_matrices, std::size_t num_rows, std::size_t num_cols)
+    : num_matrices(num_matrices),
+      num_rows(num_rows),
+      num_cols(num_cols) {
+    auto size = num_matrices * num_rows * num_cols;
+    try {
+      data = new double[size];
+    } catch (std::bad_alloc &) {
+      auto size_bytes = size * sizeof(double);
+      std::println("Memory failure! (Requested {} bytes)", size_bytes);
+      exit(EXIT_FAILURE);
+    }
+  }
+
+  tensor::tensor(const tensor &other)
+    : num_matrices(other.num_matrices),
+      num_rows(other.num_rows),
+      num_cols(other.num_cols),
+      data(other.data) {}
+
+  tensor::tensor(tensor &&other) noexcept
+    : num_matrices(other.num_matrices),
+      num_rows(other.num_rows),
+      num_cols(other.num_cols),
+      data(std::exchange(other.data, nullptr)) {}
+
+  tensor &tensor::operator=(const tensor &other) {
+    return *this = tensor(other);
+  }
+
+  tensor &tensor::operator=(tensor &&other) noexcept // move assignment
+  {
+    std::swap(data, other.data);
+    num_matrices = other.num_matrices;
+    num_cols = other.num_cols;
+    num_rows = other.num_rows;
+    return *this;
+  }
+
+  tensor::~tensor() {
+    if (data) {
+      delete[] data;
+      data = nullptr;
+    }
+  }
+
+  double &tensor::operator[](std::size_t matrix, std::size_t row, std::size_t col) {
+    return data[(num_cols * num_rows * matrix) + (num_cols * row) + (col)];
+  }
+
+  double tensor::operator[](std::size_t matrix, std::size_t row, std::size_t col) const {
+    return data[(num_cols * num_rows * matrix) + (num_cols * row) + (col)];
+  }
+
+} // namespace partdiff
diff --git a/tensor.hpp b/tensor.hpp
new file mode 100644
index 0000000..65a5992
--- /dev/null
+++ b/tensor.hpp
@@ -0,0 +1,24 @@
+#pragma once
+
+#include <cstddef>
+
+namespace partdiff {
+
+  class tensor {
+    public:
+    tensor() {};
+    tensor(std::size_t num_matrices, std::size_t num_rows, std::size_t num_cols);
+    tensor(const tensor &other);
+    tensor(tensor &&other) noexcept;
+    tensor &operator=(const tensor &other);
+    tensor &operator=(tensor &&other) noexcept;
+    ~tensor();
+    double &operator[](std::size_t matrix, std::size_t row, std::size_t col);
+    double operator[](std::size_t matrix, std::size_t row, std::size_t col) const;
+
+    private:
+    std::size_t num_matrices, num_rows, num_cols;
+    double *data = nullptr;
+  };
+
+} // namespace partdiff