diff --git a/.gitignore b/.gitignore
index 6fc13f3bfa..ffd1a8ca8d 100644
--- a/.gitignore
+++ b/.gitignore
@@ -79,6 +79,9 @@ _build
 
 *.log
 
+*.lock
+conda/recipe/meta.yaml
+conda/recipe
 .vs/
 *.key
 
@@ -89,3 +92,6 @@ uv.lock
 *.bck
 *log
 *.claude
+.pixi/
+conda/recipe/meta.yaml
+conda/recipe/meta.yaml
diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 0f3f068c29..49a5b60dca 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -2,6 +2,11 @@
 #
 # SPDX-License-Identifier: Unlicense
 
+stages:
+  - build
+  - test
+  - deploy
+
 # ============================================================
 #  GPU-specific images and runner configuration
 # ============================================================
@@ -76,15 +81,38 @@
 #  Job templates
 # ============================================================
 .build-and-test:
+  variables:
+    CMAKE_BUILD_PARALLEL_LEVEL: "3"
   script:
     - ./ci/lrz-gitlab/build-and-test.sh
 
 .benchmark:
   extends:
     - .setup-for-benchmarking
+  variables:
+    CMAKE_BUILD_PARALLEL_LEVEL: "3"
   script:
     - ./ci/lrz-gitlab/benchmark.sh
 
+.build-python-wheel:
+  variables:
+    CMAKE_BUILD_PARALLEL_LEVEL: "3"
+  before_script:
+    - python3 -m pip install --user --break-system-packages build
+  script:
+    - python3 -m build --wheel
+    - python3 ci/smoke_test_wheel.py
+  artifacts:
+    paths:
+      - dist/*.whl
+
+.pixi-python-package:
+  image: ghcr.io/prefix-dev/pixi:latest
+  variables:
+    CMAKE_BUILD_PARALLEL_LEVEL: "1"
+  before_script:
+    - pixi install --locked
+
 # ============================================================
 #  NeoN - testing jobs
 # ============================================================
@@ -126,3 +154,42 @@ benchmark-neon-intel:
     - .gpu-intel
     - .rules-neon-benchmarking
     - .benchmark
+
+# ============================================================
+#  NeoN - Python packaging job
+# ============================================================
+build-python-wheel-nvidia:
+  extends:
+    - .gpu-nvidia
+    - .rules-neon-testing
+    - .build-python-wheel
+
+build-python-wheel-pixi:
+  extends:
+    - .pixi-python-package
+  stage: build
+  script:
+    - pixi run -e build smoke-test
+  artifacts:
+    paths:
+      - dist/*.whl
+    expire_in: 1 week
+  rules:
+    - if: '$CI_COMMIT_TAG'
+    - if: '$CI_PIPELINE_SOURCE == "merge_request_event"'
+    - if: '$CI_COMMIT_BRANCH'
+
+publish-python-wheel-pypi:
+  extends:
+    - .pixi-python-package
+  stage: deploy
+  needs:
+    - job: build-python-wheel-pixi
+      artifacts: true
+  id_tokens:
+    PYPI_ID_TOKEN:
+      aud: pypi
+  script:
+    - pixi run -e release twine upload dist/*
+  rules:
+    - if: '$CI_COMMIT_TAG =~ /^v?[0-9]+\.[0-9]+\.[0-9]+$/'
diff --git a/README.md b/README.md
index 7c692e2922..89a7374889 100644
--- a/README.md
+++ b/README.md
@@ -63,6 +63,24 @@ We provide a set of unit tests which can be executed via ctest or
 
     cmake --build . --target test
 
+### Pixi Workspace
+
+For reproducible local development environments, NeoN can also be driven via
+`pixi.toml`. Pixi is used here as the environment and task runner, while the
+actual compilation and wheel creation still go through CMake and
+`python -m build`.
+
+Typical commands are:
+
+    pixi install                 # creates pixi.lock
+    pixi run -e build build      # configure + build with -j3
+    pixi run -e py312 test       # run ctest in the Python 3.12 environment
+    pixi run -e build wheel      # build a wheel without build isolation
+    pixi run -e build smoke-test # install the built wheel in an isolated venv and run test_neon_package.py
+
+The generated `pixi.lock` should be committed so CI uses pinned tooling across
+platforms.
+
 
 ## Integration with other CFD Frameworks
 
diff --git a/conda/recipe/README.md b/conda/recipe/README.md
new file mode 100644
index 0000000000..584cb14ade
--- /dev/null
+++ b/conda/recipe/README.md
@@ -0,0 +1,43 @@
+# Conda recipe for NeoN
+
+This recipe packages the NeoN Python bindings and ships the compiled extension
+(`neon/_neon*.so` on Linux/macOS, `neon/_neon*.pyd` on Windows).
+
+## Build locally
+
+```bash
+conda build conda/recipe
+```
+
+The recipe defaults to `CMAKE_BUILD_PARALLEL_LEVEL=2` to reduce memory pressure.
+Override it when needed:
+
+```bash
+CMAKE_BUILD_PARALLEL_LEVEL=1 conda build conda/recipe
+```
+
+Get the produced package path:
+
+```bash
+conda build conda/recipe --output
+```
+
+## Upload to anaconda.org
+
+```bash
+anaconda login
+anaconda upload "$(conda build conda/recipe --output)"
+```
+
+## Install and verify
+
+```bash
+conda install -c <your-channel> neon
+python -c "import neon; import neon._neon; print(neon.__version__)"
+```
+
+## Notes
+
+- The package build uses `pip install . --no-build-isolation --no-deps` with `scikit-build-core`.
+- Build requirements use Conda compilers to produce relocatable binaries.
+- Additional system/compiler requirements may be needed depending on your platform and accelerator setup.
diff --git a/conda/recipe/meta.yaml b/conda/recipe/meta.yaml
new file mode 100644
index 0000000000..a95871d838
--- /dev/null
+++ b/conda/recipe/meta.yaml
@@ -0,0 +1,49 @@
+{% set name = "exasim-neon" %}
+{% set version = "0.1.0" %}
+
+package:
+  name: {{ name|lower }}
+  version: {{ version }}
+
+source:
+  path: ../../
+
+build:
+  number: 0
+  script:
+    - export CMAKE_BUILD_PARALLEL_LEVEL=${CMAKE_BUILD_PARALLEL_LEVEL:-1}
+    - {{ PYTHON }} -m pip install . -vv --no-build-isolation --no-deps
+
+requirements:
+  build:
+    - {{ compiler('cxx') }}
+    - cmake >=3.22
+    - ninja
+  host:
+    - python
+    - pip
+    - scikit-build-core >=0.11.0
+    - nanobind >=2.9.0
+    - nlohmann_json
+    - setuptools
+  run:
+    - python
+    - numpy >=1.19.0
+
+test:
+  requires:
+    - pip
+  imports:
+    - neon
+  commands:
+    - python -c "import neon; print(neon.__version__, neon.__has_serial__)"
+    - python -c "import neon._neon as ext; print('loaded', ext.__name__)"
+
+about:
+  home: https://github.com/exasim-project/NeoN
+  summary: Python bindings for the NeoN CFD framework
+  license: MIT
+
+extra:
+  recipe-maintainers:
+    - exasim-project
diff --git a/doc/notebooks/building_blocks.ipynb b/doc/notebooks/building_blocks.ipynb
new file mode 100644
index 0000000000..574000718a
--- /dev/null
+++ b/doc/notebooks/building_blocks.ipynb
@@ -0,0 +1,236 @@
+{
+  "cells": [
+    {
+      "cell_type": "markdown",
+      "id": "19cbe4b4",
+      "metadata": {},
+      "source": [
+        "# NeoN Python demo: vectors, interpolation, and Courant number\n",
+        "\n",
+        "This second notebook demonstrates different Python capabilities of the installed `neon` package:\n",
+        "\n",
+        "1. numpy interop with NeoN vectors,\n",
+        "2. token-driven interpolation in Python,\n",
+        "3. Courant number computation from mesh + face flux.\n"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "id": "62f20ddf",
+      "metadata": {},
+      "source": [
+        "## 1) Imports and init"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": 3,
+      "id": "796d7c7f",
+      "metadata": {},
+      "outputs": [
+        {
+          "name": "stdout",
+          "output_type": "stream",
+          "text": [
+            "executor: SerialExecutor\n",
+            "cells: 100\n"
+          ]
+        }
+      ],
+      "source": [
+        "import numpy as np\n",
+        "import neon\n",
+        "\n",
+        "if not globals().get(\"_neon_initialized\", False):\n",
+        "    neon.initialize()\n",
+        "    _neon_initialized = True\n",
+        "\n",
+        "exec = neon.SerialExecutor()\n",
+        "mesh2d = neon.create_2d_uniform_mesh(exec, 10, 10)\n",
+        "\n",
+        "print(\"executor:\", exec.name())\n",
+        "print(\"cells:\", mesh2d.n_cells())"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "id": "053ed9fe",
+      "metadata": {},
+      "source": [
+        "## 2) Capability: numpy view over NeoN vectors\n",
+        "\n",
+        "`np.asarray(neon_vector)` gives a numpy view (CPU executors), which is handy for quick analysis and plotting."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": 8,
+      "id": "953fc8f6",
+      "metadata": {},
+      "outputs": [
+        {
+          "name": "stdout",
+          "output_type": "stream",
+          "text": [
+            "initial: [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]\n",
+            "after numpy edit: [1.  1.  1.  4.5 4.5 4.5 4.5 1.  1.  1. ]\n",
+            "all values still 1.0? False\n"
+          ]
+        }
+      ],
+      "source": [
+        "v = neon.ScalarVector(exec, 10, 1.0)\n",
+        "a = np.asarray(v)\n",
+        "\n",
+        "print(\"initial:\", a)\n",
+        "\n",
+        "a[3:7] = 4.5\n",
+        "print(\"after numpy edit:\", a)\n",
+        "print(\"all values still 1.0?\", neon.equal(v, 1.0))"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "id": "5d2e42b2",
+      "metadata": {},
+      "source": [
+        "## 3) Capability: vector field interpolation selected by runtime tokens\n",
+        "\n",
+        "We use `TokenList` to pick interpolation scheme at runtime."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": 5,
+      "id": "c1d2b5c9",
+      "metadata": {},
+      "outputs": [
+        {
+          "name": "stdout",
+          "output_type": "stream",
+          "text": [
+            "surface field shape: (220, 3)\n",
+            "first 3 face vectors:\n",
+            " [[1. 2. 3.]\n",
+            " [1. 2. 3.]\n",
+            " [1. 2. 3.]]\n"
+          ]
+        }
+      ],
+      "source": [
+        "U = neon.VectorVolumeField(exec, \"U\", mesh2d)\n",
+        "neon.fill(U.internal_vector(), neon.Vec3(1.0, 2.0, 3.0))\n",
+        "\n",
+        "tokens = neon.TokenList()\n",
+        "tokens.insert_string(\"linear\")\n",
+        "\n",
+        "interp = neon.SurfaceInterpolationVec3(exec, mesh2d, tokens)\n",
+        "U_face = interp.interpolate(U)\n",
+        "U_face_np = np.asarray(U_face.internal_vector())\n",
+        "\n",
+        "print(\"surface field shape:\", U_face_np.shape)\n",
+        "print(\"first 3 face vectors:\\n\", U_face_np[:3])\n"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "id": "7678ea63",
+      "metadata": {},
+      "source": [
+        "## 4) Capability: Courant number from mesh + face flux\n",
+        "\n",
+        "`compute_co_num` returns `(max_Co, mean_Co)`."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": 6,
+      "id": "3e1ca1d9",
+      "metadata": {},
+      "outputs": [
+        {
+          "ename": "NameError",
+          "evalue": "name 'mesh' is not defined",
+          "output_type": "error",
+          "traceback": [
+            "\u001b[31m---------------------------------------------------------------------------\u001b[39m",
+            "\u001b[31mNameError\u001b[39m                                 Traceback (most recent call last)",
+            "\u001b[36mCell\u001b[39m\u001b[36m \u001b[39m\u001b[32mIn[6]\u001b[39m\u001b[32m, line 1\u001b[39m\n\u001b[32m----> \u001b[39m\u001b[32m1\u001b[39m face_count = mesh2d.n_internal_faces() + \u001b[43mmesh\u001b[49m.n_boundary_faces()\n\u001b[32m      2\u001b[39m face_flux = neon.ScalarVector(exec, face_count)\n\u001b[32m      3\u001b[39m neon.fill(face_flux, \u001b[32m1.0\u001b[39m)\n",
+            "\u001b[31mNameError\u001b[39m: name 'mesh' is not defined"
+          ]
+        }
+      ],
+      "source": [
+        "face_count = mesh2d.n_internal_faces() + mesh.n_boundary_faces()\n",
+        "face_flux = neon.ScalarVector(exec, face_count)\n",
+        "neon.fill(face_flux, 1.0)\n",
+        "\n",
+        "dt = 0.01\n",
+        "max_co, mean_co = neon.compute_co_num(mesh, face_flux, dt)\n",
+        "\n",
+        "print(\"max Co:\", max_co)\n",
+        "print(\"mean Co:\", mean_co)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "id": "5a6fba33",
+      "metadata": {},
+      "source": [
+        "## 5) Optional finalize"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": 12,
+      "id": "53dc543a",
+      "metadata": {},
+      "outputs": [
+        {
+          "name": "stdout",
+          "output_type": "stream",
+          "text": [
+            "NeoN finalized\n"
+          ]
+        }
+      ],
+      "source": [
+        "if globals().get(\"_neon_initialized\", False):\n",
+        "    neon.finalize()\n",
+        "    _neon_initialized = False\n",
+        "    print(\"NeoN finalized\")\n",
+        "else:\n",
+        "    print(\"NeoN already finalized\")"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "id": "e14ea8eb-9c86-4356-b80c-06ca839113c4",
+      "metadata": {},
+      "outputs": [],
+      "source": []
+    }
+  ],
+  "metadata": {
+    "kernelspec": {
+      "display_name": "Python 3 (ipykernel)",
+      "language": "python",
+      "name": "python3"
+    },
+    "language_info": {
+      "codemirror_mode": {
+        "name": "ipython",
+        "version": 3
+      },
+      "file_extension": ".py",
+      "mimetype": "text/x-python",
+      "name": "python",
+      "nbconvert_exporter": "python",
+      "pygments_lexer": "ipython3",
+      "version": "3.13.1"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 5
+}
diff --git a/doc/notebooks/expression_example.ipynb b/doc/notebooks/expression_example.ipynb
new file mode 100644
index 0000000000..dc2010d681
--- /dev/null
+++ b/doc/notebooks/expression_example.ipynb
@@ -0,0 +1,450 @@
+{
+  "cells": [
+    {
+      "cell_type": "markdown",
+      "id": "9880b8a3",
+      "metadata": {},
+      "source": [
+        "# NeoN Python coverage: benchmark `laplacian(gamma, U) - div(phi, U)`\n",
+        "\n",
+        "This notebook now tracks every major element in the C++ benchmark snippet:\n",
+        "1. executor sweep and size sweep metadata,\n",
+        "2. `nCells = 10` mesh and `U`, `phi`, `gamma` setup,\n",
+        "3. unoptimized and fused (optimized) expression paths,\n",
+        "4. runtime scheme tokens (`Gauss`, `linear`, `uncorrected`, `upwind`),\n",
+        "5. benchmark-style loops over executors and sizes.\n",
+        "\n",
+        "> Current Python bindings expose operator construction but not `Expression.read(...)`, `Expression.assemble(...)`, or `dsl::optimize(...)`. This notebook explicitly demonstrates parity where available and marks the non-exposed steps."
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "id": "e902f751",
+      "metadata": {},
+      "source": [
+        "## 1) Imports and runtime init"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": 1,
+      "id": "3c2cebd1",
+      "metadata": {},
+      "outputs": [
+        {
+          "name": "stdout",
+          "output_type": "stream",
+          "text": [
+            "cwd: /Volumes/Data/Code/NeoN/doc/notebooks\n",
+            "repo_root: /Volumes/Data/Code/NeoN\n",
+            "local_bindings: /Volumes/Data/Code/NeoN/build/develop/bindings\n",
+            "neon module: /Volumes/Data/Code/NeoN/build/develop/bindings/neon/__init__.py\n",
+            "executor: SerialExecutor\n",
+            "cells: 8\n"
+          ]
+        }
+      ],
+      "source": [
+        "import os\n",
+        "import sys\n",
+        "from pathlib import Path\n",
+        "import numpy as np\n",
+        "\n",
+        "def _find_repo_root(start: Path) -> Path | None:\n",
+        "    for candidate in [start, *start.parents]:\n",
+        "        if (candidate / \"CMakeLists.txt\").exists() and (candidate / \"doc\" / \"notebooks\").exists():\n",
+        "            return candidate\n",
+        "    return None\n",
+        "\n",
+        "cwd = Path(os.getcwd()).resolve()\n",
+        "repo_root = _find_repo_root(cwd)\n",
+        "if repo_root is None:\n",
+        "    fallback = Path(\"/home/andrei2/Downloads/NeoN\").resolve()\n",
+        "    if (fallback / \"CMakeLists.txt\").exists():\n",
+        "        repo_root = fallback\n",
+        "    else:\n",
+        "        raise RuntimeError(\"Could not locate NeoN repo root from notebook working directory\")\n",
+        "\n",
+        "local_bindings = repo_root / \"build\" / \"develop\" / \"bindings\"\n",
+        "if not local_bindings.exists():\n",
+        "    raise RuntimeError(f\"Local bindings folder not found: {local_bindings}\")\n",
+        "\n",
+        "if str(local_bindings) in sys.path:\n",
+        "    sys.path.remove(str(local_bindings))\n",
+        "sys.path.insert(0, str(local_bindings))\n",
+        "\n",
+        "for mod_name in [m for m in list(sys.modules) if m == \"neon\" or m.startswith(\"neon.\")]:\n",
+        "    del sys.modules[mod_name]\n",
+        "\n",
+        "import neon\n",
+        "from neon import imp, exp\n",
+        "\n",
+        "if not globals().get(\"_neon_initialized\", False):\n",
+        "    neon.initialize()\n",
+        "    _neon_initialized = True\n",
+        "\n",
+        "exec = neon.SerialExecutor()\n",
+        "mesh = neon.create_1d_uniform_mesh(exec, 8)\n",
+        "\n",
+        "print(\"cwd:\", cwd)\n",
+        "print(\"repo_root:\", repo_root)\n",
+        "print(\"local_bindings:\", local_bindings)\n",
+        "print(\"neon module:\", neon.__file__)\n",
+        "print(\"executor:\", exec.name())\n",
+        "print(\"cells:\", mesh.n_cells())"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "id": "3f27f181",
+      "metadata": {},
+      "source": [
+        "## 2) Benchmark metadata: sizes, epsilon, and available executors"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": 2,
+      "id": "a9de755f",
+      "metadata": {},
+      "outputs": [
+        {
+          "name": "stdout",
+          "output_type": "stream",
+          "text": [
+            "epsilon: 1e-32\n",
+            "sizes: [65536, 131072, 262144, 524288, 1048576]\n",
+            "executors: ['SerialExecutor', 'CPUExecutor']\n"
+          ]
+        }
+      ],
+      "source": [
+        "epsilon = 1e-32\n",
+        "sizes = [1 << 16, 1 << 17, 1 << 18, 1 << 19, 1 << 20]\n",
+        "\n",
+        "executors = [(\"SerialExecutor\", neon.SerialExecutor()), (\"CPUExecutor\", neon.CPUExecutor())]\n",
+        "if neon.gpu_available():\n",
+        "    executors.append((\"GPUExecutor\", neon.GPUExecutor()))\n",
+        "\n",
+        "print(\"epsilon:\", epsilon)\n",
+        "print(\"sizes:\", sizes)\n",
+        "print(\"executors:\", [name for name, _ in executors])"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "id": "3e3767a5",
+      "metadata": {},
+      "source": [
+        "## 3) Mirror C++ field setup (`nCells = 10`, `nFaces = 9`)"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": 4,
+      "id": "442ddbd9",
+      "metadata": {},
+      "outputs": [
+        {
+          "name": "stdout",
+          "output_type": "stream",
+          "text": [
+            "executor: SerialExecutor\n",
+            "nCells: 10\n",
+            "nFaces (mesh): 11 (benchmark constant is 9)\n",
+            "volume BC count: 2 surface BC count: 2\n"
+          ]
+        }
+      ],
+      "source": [
+        "exec_name, exec = executors[0]\n",
+        "n_cells = 10\n",
+        "mesh = neon.create_1d_uniform_mesh(exec, n_cells)\n",
+        "\n",
+        "vol_bcs = neon.create_calculated_volume_bcs_scalar(mesh)\n",
+        "surf_bcs = neon.create_calculated_surface_bcs_scalar(mesh)\n",
+        "\n",
+        "U = neon.ScalarVolumeField(exec, \"U\", mesh)\n",
+        "phi = neon.ScalarSurfaceField(exec, \"phi\", mesh)\n",
+        "gamma = neon.ScalarSurfaceField(exec, \"gamma\", mesh)\n",
+        "\n",
+        "neon.fill(U.internal_vector(), 2.0)\n",
+        "neon.fill(phi.internal_vector(), 1.0)\n",
+        "neon.fill(gamma.internal_vector(), 2.0)\n",
+        "\n",
+        "n_faces = mesh.n_internal_faces() + mesh.n_boundary_faces()\n",
+        "print(\"executor:\", exec_name)\n",
+        "print(\"nCells:\", mesh.n_cells())\n",
+        "print(\"nFaces (mesh):\", n_faces, \"(benchmark constant is 9)\")\n",
+        "print(\"volume BC count:\", len(vol_bcs), \"surface BC count:\", len(surf_bcs))\n"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "id": "0fa01bdb",
+      "metadata": {},
+      "source": [
+        "## 4) Unoptimized expression path (`laplacian - div`)"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": 5,
+      "id": "d2a9d153",
+      "metadata": {},
+      "outputs": [
+        {
+          "name": "stdout",
+          "output_type": "stream",
+          "text": [
+            "laplacian operator: LaplacianOperator\n",
+            "divergence operator: DivOperator\n",
+            "unoptimized expression terms: 2\n"
+          ]
+        }
+      ],
+      "source": [
+        "lap_op = imp.laplacian(gamma, U)\n",
+        "div_op = imp.div(phi, U)\n",
+        "expr_unoptimized = lap_op - div_op\n",
+        "\n",
+        "print(\"laplacian operator:\", lap_op.get_name())\n",
+        "print(\"divergence operator:\", div_op.get_name())\n",
+        "print(\"unoptimized expression terms:\", expr_unoptimized.size())"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "id": "02ed89bc",
+      "metadata": {},
+      "source": [
+        "## 5) Fused path mapping (`dsl::optimize(expr)` in C++)"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": 6,
+      "id": "b0a5bea7",
+      "metadata": {},
+      "outputs": [
+        {
+          "name": "stdout",
+          "output_type": "stream",
+          "text": [
+            "Python exposes neon.optimize: False\n",
+            "using expr_unoptimized as fused placeholder (binding not exposed yet)\n"
+          ]
+        }
+      ],
+      "source": [
+        "expr_fused = expr_unoptimized\n",
+        "optimize_exposed = hasattr(neon, \"optimize\")\n",
+        "\n",
+        "print(\"Python exposes neon.optimize:\", optimize_exposed)\n",
+        "if optimize_exposed:\n",
+        "    expr_fused = neon.optimize(expr_unoptimized)\n",
+        "    print(\"fused expression terms:\", expr_fused.size())\n",
+        "else:\n",
+        "    print(\"using expr_unoptimized as fused placeholder (binding not exposed yet)\")\n"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "id": "036f77f1",
+      "metadata": {},
+      "source": [
+        "## 6) Scheme-token coverage (`Gauss`, `linear`, `uncorrected`, `upwind`)"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": 7,
+      "id": "fd334d81",
+      "metadata": {},
+      "outputs": [
+        {
+          "name": "stdout",
+          "output_type": "stream",
+          "text": [
+            "laplacian(gamma,U) tokens: ['Gauss', 'linear', 'uncorrected']\n",
+            "div(phi,U) tokens: ['Gauss', 'upwind']\n",
+            "read(Dictionary) coverage note: nested Dictionary->TokenList scheme trees are not exposed in Python yet\n"
+          ]
+        }
+      ],
+      "source": [
+        "lap_tokens = neon.TokenList()\n",
+        "lap_tokens.insert_string(\"Gauss\")\n",
+        "lap_tokens.insert_string(\"linear\")\n",
+        "lap_tokens.insert_string(\"uncorrected\")\n",
+        "\n",
+        "div_tokens = neon.TokenList()\n",
+        "div_tokens.insert_string(\"Gauss\")\n",
+        "div_tokens.insert_string(\"upwind\")\n",
+        "\n",
+        "print(\"laplacian(gamma,U) tokens:\", [lap_tokens.get_string(i) for i in range(lap_tokens.size())])\n",
+        "print(\"div(phi,U) tokens:\", [div_tokens.get_string(i) for i in range(div_tokens.size())])\n",
+        "\n",
+        "print(\"read(Dictionary) coverage note: nested Dictionary->TokenList scheme trees are not exposed in Python yet\")\n"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "id": "86e31369",
+      "metadata": {},
+      "source": [
+        "## 7) Benchmark-style loop coverage (executors x sizes)"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": 8,
+      "id": "4a4b986d",
+      "metadata": {},
+      "outputs": [
+        {
+          "name": "stdout",
+          "output_type": "stream",
+          "text": [
+            "rows produced: 10\n",
+            "first 5 rows:\n",
+            "('SerialExecutor', 65536, 2, 2)\n",
+            "('SerialExecutor', 131072, 2, 2)\n",
+            "('SerialExecutor', 262144, 2, 2)\n",
+            "('SerialExecutor', 524288, 2, 2)\n",
+            "('SerialExecutor', 1048576, 2, 2)\n",
+            "assemble coverage note: Expression.assemble(...) is now exposed (with DB-registered fields)\n"
+          ]
+        }
+      ],
+      "source": [
+        "rows = []\n",
+        "for exec_name, exec in executors:\n",
+        "    mesh = neon.create_1d_uniform_mesh(exec, 10)\n",
+        "    U = neon.ScalarVolumeField(exec, \"U\", mesh)\n",
+        "    phi = neon.ScalarSurfaceField(exec, \"phi\", mesh)\n",
+        "    gamma = neon.ScalarSurfaceField(exec, \"gamma\", mesh)\n",
+        "    neon.fill(U.internal_vector(), 2.0)\n",
+        "    neon.fill(phi.internal_vector(), 1.0)\n",
+        "    neon.fill(gamma.internal_vector(), 2.0)\n",
+        "\n",
+        "    expr_unopt = imp.laplacian(gamma, U) - imp.div(phi, U)\n",
+        "    expr_fused = neon.optimize(expr_unopt) if hasattr(neon, \"optimize\") else expr_unopt\n",
+        "\n",
+        "    for size in sizes:\n",
+        "        rows.append((exec_name, size, expr_unopt.size(), expr_fused.size()))\n",
+        "\n",
+        "print(\"rows produced:\", len(rows))\n",
+        "print(\"first 5 rows:\")\n",
+        "for row in rows[:5]:\n",
+        "    print(row)\n",
+        "\n",
+        "print(\"assemble coverage note: Expression.assemble(...) is now exposed (with DB-registered fields)\")\n"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "id": "e453934b",
+      "metadata": {},
+      "source": [
+        "## 8) New: implicit assembly from Python (`register_volume_field` + `assemble`)"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": 9,
+      "id": "96dd8dba",
+      "metadata": {},
+      "outputs": [
+        {
+          "name": "stdout",
+          "output_type": "stream",
+          "text": [
+            "expression terms: 2\n",
+            "sparsity rows: 10\n",
+            "sparsity nnz: 28\n",
+            "linear system type: LinearSystemScalar\n"
+          ]
+        }
+      ],
+      "source": [
+        "exec = neon.SerialExecutor()\n",
+        "mesh = neon.create_1d_uniform_mesh(exec, 10)\n",
+        "\n",
+        "phi = neon.ScalarVolumeField(exec, \"phi\", mesh)\n",
+        "neon.fill(phi.internal_vector(), 1.0)\n",
+        "\n",
+        "db = neon.Database()\n",
+        "phi_reg = neon.register_volume_field(db, \"fields\", phi)\n",
+        "\n",
+        "eq = imp.source(phi_reg, phi_reg) + imp.source(phi_reg, phi_reg)\n",
+        "sparsity, linear_system = eq.assemble(mesh, 0.0, 1.0)\n",
+        "\n",
+        "print(\"expression terms:\", eq.size())\n",
+        "print(\"sparsity rows:\", sparsity.rows())\n",
+        "print(\"sparsity nnz:\", sparsity.nnz())\n",
+        "print(\"linear system type:\", type(linear_system).__name__)\n"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "id": "c716a35c",
+      "metadata": {},
+      "source": [
+        "## 9) Optional finalize"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": 10,
+      "id": "6eba9dae",
+      "metadata": {},
+      "outputs": [
+        {
+          "name": "stdout",
+          "output_type": "stream",
+          "text": [
+            "NeoN finalized\n"
+          ]
+        }
+      ],
+      "source": [
+        "if globals().get(\"_neon_initialized\", False):\n",
+        "    neon.finalize()\n",
+        "    _neon_initialized = False\n",
+        "    print(\"NeoN finalized\")\n",
+        "else:\n",
+        "    print(\"NeoN already finalized\")\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "id": "c228f791-ba45-4f6a-a3c3-dde9b3a2f368",
+      "metadata": {},
+      "outputs": [],
+      "source": []
+    }
+  ],
+  "metadata": {
+    "kernelspec": {
+      "display_name": "Python 3 (ipykernel)",
+      "language": "python",
+      "name": "python3"
+    },
+    "language_info": {
+      "codemirror_mode": {
+        "name": "ipython",
+        "version": 3
+      },
+      "file_extension": ".py",
+      "mimetype": "text/x-python",
+      "name": "python",
+      "nbconvert_exporter": "python",
+      "pygments_lexer": "ipython3",
+      "version": "3.13.1"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 5
+}
diff --git a/pyproject.toml b/pyproject.toml
index 53c3cff4c6..25bff43005 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -3,7 +3,7 @@
 # SPDX-License-Identifier: Unlicense
 
 [project]
-name = "neon_pde"
+name = "exasim-neon"
 version = "0.1.0"
 authors = [
     {name = "NeoN authors", email = "neon@example.com"}
@@ -40,7 +40,7 @@ test = [
     "pytest>=6.0",
 ]
 all = [
-    "neon_pde[dev,test]",
+    "exasim-neon[dev,test]",
 ]
 
 [build-system]
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 0000000000..84933bd9e7
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,2 @@
+numpy>=1.19.0
+pytest>=6.0
diff --git a/src/bindings/CMakeLists.txt b/src/bindings/CMakeLists.txt
index fe69b88b34..d58bd23d1f 100644
--- a/src/bindings/CMakeLists.txt
+++ b/src/bindings/CMakeLists.txt
@@ -22,8 +22,8 @@ if(DEFINED SKBUILD)
   # scikit-build-core manages installation paths
   set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR})
 else()
-  # Standalone build - convention is to put bindings .so file in a project name folder
-  set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bindings/neon)
+  # Standalone build - always place bindings in a fixed path for tests and local usage
+  set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}/build/develop/bindings/neon)
 endif()
 
 nanobind_add_module(
@@ -105,9 +105,9 @@ if(DEFINED SKBUILD)
   install(FILES ${CMAKE_CURRENT_SOURCE_DIR}/../neon/py.typed DESTINATION neon)
 else()
   # Standalone installation
-  install(TARGETS _neon LIBRARY DESTINATION ${CMAKE_INSTALL_LIBDIR}/python/neon)
-  install(FILES ${CMAKE_LIBRARY_OUTPUT_DIRECTORY}/__init__.py
-          DESTINATION ${CMAKE_INSTALL_LIBDIR}/python/neon)
+  install(TARGETS _neon LIBRARY DESTINATION ${Python_SITELIB}/neon)
+  install(FILES ${CMAKE_LIBRARY_OUTPUT_DIRECTORY}/__init__.py DESTINATION ${Python_SITELIB}/neon)
   install(FILES ${CMAKE_CURRENT_SOURCE_DIR}/../neon/py.typed
           DESTINATION ${CMAKE_INSTALL_LIBDIR}/python/neon)
+
 endif()
diff --git a/src/bindings/dsl.cpp b/src/bindings/dsl.cpp
index 31648387b0..e28ce6a222 100644
--- a/src/bindings/dsl.cpp
+++ b/src/bindings/dsl.cpp
@@ -4,6 +4,7 @@
 
 #include <nanobind/nanobind.h>
 #include <nanobind/stl/string.h>
+#include <nanobind/stl/tuple.h>
 #include <nanobind/stl/vector.h>
 #include <nanobind/stl/variant.h>
 #include <nanobind/operators.h>
@@ -135,6 +136,17 @@ void declare_dsl_components(nb::module_& m, const std::string& suffix)
         .def(
             "__sub__", [](Expr lhs, const TemporalOp& rhs) { return lhs - rhs; }, nb::is_operator()
         )
+        .def(
+            "assemble",
+            [](const Expr& expr, const UnstructuredMesh& mesh, scalar t, scalar dt)
+            {
+                auto [sparsity, linear_system] = expr.assemble(mesh, t, dt);
+                return std::make_tuple(*sparsity, std::move(linear_system));
+            },
+            "mesh"_a,
+            "t"_a,
+            "dt"_a
+        )
         .def("size", &Expr::size);
 }
 
diff --git a/src/bindings/unstructuredMesh.cpp b/src/bindings/unstructuredMesh.cpp
index d7e2a10926..a4746570fd 100644
--- a/src/bindings/unstructuredMesh.cpp
+++ b/src/bindings/unstructuredMesh.cpp
@@ -197,6 +197,42 @@ void registerUnstructuredMesh(nb::module_& m)
         "Each cell has a left and right face. Useful for 1D simulations\n"
         "and testing finite volume schemes."
     );
+    m.def(
+        "create_2d_uniform_mesh",
+        &NeoN::create2DUniformMesh,
+        "exec"_a,
+        "nx"_a,
+        "ny"_a,
+        "Lx"_a = 1.0,
+        "Ly"_a = 1.0,
+        "Create a uniform 2D mesh on [0,Lx] x [0,Ly].\n\n"
+        "Args:\n"
+        "    exec: Executor for parallel operations\n"
+        "    nx: Number of cells in the x-direction\n"
+        "    ny: Number of cells in the y-direction\n"
+        "    Lx: Length in the x-direction (Default: 1.0)\n"
+        "    Ly: Length in the y-direction (Default: 1.0)"
+    );
+    m.def(
+        "create_3d_uniform_mesh",
+        &NeoN::create3DUniformMesh,
+        "exec"_a,
+        "nx"_a,
+        "ny"_a,
+        "nz"_a,
+        "Lx"_a = 1.0,
+        "Ly"_a = 1.0,
+        "Lz"_a = 1.0,
+        "Create a uniform 3D mesh on [0,Lx] x [0,Ly] x [0,Lz].\n\n"
+        "Args:\n"
+        "    exec: Executor for parallel operations\n"
+        "    nx: Number of cells in the x-direction\n"
+        "    ny: Number of cells in the y-direction\n"
+        "    nz: Number of cells in the z-direction\n"
+        "    Lx: Length in the x-direction (Default: 1.0)\n"
+        "    Ly: Length in the y-direction (Default: 1.0)\n"
+        "    Lz: Length in the z-direction (Default: 1.0)"
+    );
 }
 
 } // namespace NeoN::bindings
diff --git a/src/bindings/volumeField.cpp b/src/bindings/volumeField.cpp
index 1b716dd529..8e332e2012 100644
--- a/src/bindings/volumeField.cpp
+++ b/src/bindings/volumeField.cpp
@@ -279,6 +279,49 @@ void registerVolumeField(nb::module_& m)
         "Rotate old-time Vec3 volume field (φ^n → φ^{n-1}) — field must be registered in "
         "VectorCollection"
     );
+
+    m.def(
+        "register_volume_field",
+        [](NeoN::Database& db,
+           const std::string& collection_name,
+           fvcc::VolumeField<NeoN::scalar>& field) -> fvcc::VolumeField<NeoN::scalar>&
+        {
+            auto& collection = fvcc::VectorCollection::instance(db, collection_name);
+            return collection.registerVector<fvcc::VolumeField<NeoN::scalar>>(
+                fvcc::CreateFromExistingVector<fvcc::VolumeField<NeoN::scalar>> {
+                    .name = field.name,
+                    .field = field,
+                }
+            );
+        },
+        "db"_a,
+        "collection_name"_a,
+        "field"_a,
+        nb::rv_policy::reference,
+        "Register a scalar volume field in a database collection and return the registered "
+        "field"
+    );
+
+    m.def(
+        "register_volume_field",
+        [](NeoN::Database& db,
+           const std::string& collection_name,
+           fvcc::VolumeField<NeoN::Vec3>& field) -> fvcc::VolumeField<NeoN::Vec3>&
+        {
+            auto& collection = fvcc::VectorCollection::instance(db, collection_name);
+            return collection.registerVector<fvcc::VolumeField<NeoN::Vec3>>(
+                fvcc::CreateFromExistingVector<fvcc::VolumeField<NeoN::Vec3>> {
+                    .name = field.name,
+                    .field = field,
+                }
+            );
+        },
+        "db"_a,
+        "collection_name"_a,
+        "field"_a,
+        nb::rv_policy::reference,
+        "Register a Vec3 volume field in a database collection and return the registered field"
+    );
 }
 
 } // namespace NeoN::bindings
diff --git a/test/bindings/CMakeLists.txt b/test/bindings/CMakeLists.txt
index 3af00ac51b..a9fd4f54da 100644
--- a/test/bindings/CMakeLists.txt
+++ b/test/bindings/CMakeLists.txt
@@ -2,6 +2,19 @@
 #
 # SPDX-License-Identifier: Unlicense
 
+find_package(Python COMPONENTS Interpreter Development)
+set(PYTHON_EXECUTABLE ${Python_EXECUTABLE})
+
+# * Function to generate bindings unit test.
+function(neon_bindings_test TEST)
+  # Correct variable according to https://cmake.org/cmake/help/latest/module/FindPython.html
+  # Printing on test skip: https://stackoverflow.com/a/64776498
+  add_test(
+    NAME neon_bindings_test_${TEST}
+    COMMAND ${Python_EXECUTABLE} -m pytest -vv -rs ${CMAKE_CURRENT_SOURCE_DIR}/test_${TEST}.py
+    WORKING_DIRECTORY ${CMAKE_BINARY_DIR}/bindings/)
+endfunction()
+
 if(NeoN_BUILD_PYTHON_BINDINGS)
   add_test(
     NAME neon_bindings_tests
diff --git a/test/bindings/conftest.py b/test/bindings/conftest.py
index 229a6d63f8..50d165d10d 100644
--- a/test/bindings/conftest.py
+++ b/test/bindings/conftest.py
@@ -5,7 +5,7 @@
 """
 Test configuration for NeoN Python bindings.
 
-This module automatically locates the neon package from the build directory
+This module imports the neon package from a fixed build path
 and configures pytest fixtures for executor parameterization.
 """
 
@@ -13,11 +13,15 @@
 
 # Set up path before importing neon
 import neon
+import sys
+import types
 
+import neon
 
 @pytest.fixture(scope="session", autouse=True)
 def neon_global_session():
     """Initialize NeoN once for all test files in this session."""
+
     neon.initialize()
     yield  # This is where all tests run
     neon.finalize()
@@ -30,6 +34,9 @@ def pytest_configure(config):
 
 def get_available_executors():
     """Get list of executors available at build time."""
+    if not NEON_AVAILABLE:
+        return []
+
     executors = []
     if neon.__has_serial__:
         executors.append(("serial", neon.SerialExecutor))
@@ -56,6 +63,7 @@ def pytest_generate_tests(metafunc):
 
 
 def pytest_collection_modifyitems(config, items):
+    """Skip tests based on runtime availability."""
     """Automatically skip GPU tests if GPU is not available."""
     if hasattr(neon, '__has_gpu__') and not neon.__has_gpu__:
         skip_gpu = pytest.mark.skip(reason="GPU not available")
diff --git a/test/bindings/test_dsl.py b/test/bindings/test_dsl.py
index 8520041e54..68bb769e08 100644
--- a/test/bindings/test_dsl.py
+++ b/test/bindings/test_dsl.py
@@ -152,3 +152,21 @@ def test_dsl_vector_operators(executor):
     res_eqn = eqn1 - eqn2
     assert isinstance(res_eqn, neon.ExpressionVector)
     assert res_eqn.size() == 4
+
+
+def test_dsl_scalar_assemble(executor):
+    name, exec = executor
+    mesh = neon.create_1d_uniform_mesh(exec, 10)
+    phi = neon.ScalarVolumeField(exec, "phi", mesh)
+    db = neon.Database()
+    phi_reg = neon.register_volume_field(db, "fields", phi)
+
+    equation = imp.source(phi_reg, phi_reg) + imp.source(phi_reg, phi_reg)
+    assert hasattr(equation, "assemble")
+
+    sparsity, linear_system = equation.assemble(mesh, 0.0, 1.0)
+
+    assert isinstance(sparsity, neon.SparsityPattern)
+    assert isinstance(linear_system, neon.LinearSystemScalar)
+    assert sparsity.rows() == mesh.n_cells()
+    assert sparsity.nnz() > 0
diff --git a/test/bindings/test_scalar.py b/test/bindings/test_scalar.py
index 7c4f4826ba..23e6ca77df 100644
--- a/test/bindings/test_scalar.py
+++ b/test/bindings/test_scalar.py
@@ -23,7 +23,7 @@ def test_scalar_imports():
 def test_rootvsmall():
     """Test ROOTVSMALL constant."""
     assert neon.ROOTVSMALL > 0
-    assert neon.ROOTVSMALL < 1e-10  # Should be very small
+    assert neon.ROOTVSMALL < 1e-10
 
 
 def test_double_precision_flag():
@@ -33,14 +33,12 @@ def test_double_precision_flag():
 
 def test_scalar_one():
     """Test scalar_one function."""
-    one = neon.scalar_one()
-    assert one == 1.0
+    assert neon.scalar_one() == 1.0
 
 
 def test_scalar_zero():
     """Test scalar_zero function."""
-    zero = neon.scalar_zero()
-    assert zero == 0.0
+    assert neon.scalar_zero() == 0.0
 
 
 def test_scalar_mag_positive():
diff --git a/test/bindings/test_vector.py b/test/bindings/test_vector.py
index 2bf1efbd77..ff4f1c9b1f 100644
--- a/test/bindings/test_vector.py
+++ b/test/bindings/test_vector.py
@@ -9,6 +9,7 @@
     has_jax = False
 
 import neon
+import numpy as np
 
 def test_scalar_vector(executor):
     name, exec = executor
@@ -68,8 +69,8 @@ def test_copy_to_host():
     assert vv1.size() == vv2.size()
 
 
-def test_jax_operations(executor):
-    """Test JAX operations on NeoN ScalarVectors.
+def test_array_operations(executor):
+    """Test array operations on NeoN ScalarVectors via numpy.
 
     For CPU executors (Serial/CPU): convert directly via __array__.
     For GPU executors: copy to host first, then convert.
@@ -81,12 +82,11 @@ def test_jax_operations(executor):
     values = [1.0, 2.0, 3.0, 4.0, 5.0]
     v = neon.ScalarVector(exec, values)
 
-    # GPU vectors must be copied to host before converting to JAX
     host_v = v.copy_to_host() if name == "gpu" else v
 
-    arr = jnp.asarray(host_v)
+    arr = np.asarray(host_v)
     assert arr.shape == (5,)
-    assert float(jnp.sum(arr)) == 15.0
-    assert float(jnp.mean(arr)) == 3.0
-    assert float(jnp.max(arr)) == 5.0
-    assert float(jnp.min(arr)) == 1.0
+    assert float(np.sum(arr)) == 15.0
+    assert float(np.mean(arr)) == 3.0
+    assert float(np.max(arr)) == 5.0
+    assert float(np.min(arr)) == 1.0
diff --git a/test_neon_package.py b/test_neon_package.py
new file mode 100644
index 0000000000..ad8849513a
--- /dev/null
+++ b/test_neon_package.py
@@ -0,0 +1,47 @@
+#!/usr/bin/env python3
+
+# SPDX-FileCopyrightText: 2026 NeoN authors
+#
+# SPDX-License-Identifier: MIT
+
+"""Smoke test for the installed NeoN Python package."""
+
+from __future__ import annotations
+
+import json
+import sys
+import traceback
+
+
+def main() -> int:
+    try:
+        import neon
+        import neon._neon as ext
+    except Exception as exc:
+        print("Failed to import NeoN package:", exc)
+        traceback.print_exc()
+        return 1
+
+    info = {
+        "version": getattr(neon, "__version__", "unknown"),
+        "extension": getattr(ext, "__name__", "unknown"),
+        "has_serial": getattr(neon, "__has_serial__", None),
+        "has_cpu": getattr(neon, "__has_cpu__", None),
+        "has_gpu": getattr(neon, "__has_gpu__", None),
+    }
+    print(json.dumps(info, indent=2))
+
+    try:
+        neon.initialize()
+        neon.finalize()
+    except Exception as exc:
+        print("Runtime check failed:", exc)
+        traceback.print_exc()
+        return 2
+
+    print("NeoN package smoke test: OK")
+    return 0
+
+
+if __name__ == "__main__":
+    sys.exit(main())