diff --git a/.gitignore b/.gitignore
index b37ba8c..4afdf81 100644
--- a/.gitignore
+++ b/.gitignore
@@ -7,6 +7,7 @@
 # These are backup files generated by rustfmt
 **/*.rs.bk
 
-# IDE
+# IDE & Agents
 .idea/
 .vscode/
+.codex
diff --git a/CHANGELOG.md b/CHANGELOG.md
index d624ff0..3d4f21a 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,8 +1,15 @@
 # Changelog
 
-## v0.4.0 (2022)
+## v0.5.0 (01.05.2026)
+
+* Fix 3D spline evaluation layout handling for multivariate data.
+* Update dependencies for `ndarray` 0.17 and replace deprecated reshape calls with explicit layout-aware reshaping.
+
+
+## v0.4.0 (20.04.2022)
 * Update to ndarray 0.15 and sprs 0.11
 
+
 ## v0.3.0 (22.03.2020)
 
 * Add `Real` pub trait. Trait `Real` is only implemented for `f32` and `f64`, 
diff --git a/Cargo.lock b/Cargo.lock
index 403f0ca..5dfb907 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -1,6 +1,6 @@
 # This file is automatically @generated by Cargo.
 # It is not intended for manual editing.
-version = 3
+version = 4
 
 [[package]]
 name = "alga"
@@ -70,7 +70,7 @@ checksum = "22ec99545bb0ed0ea7bb9b8e1e9122ea386ff8a48c0922e43f36d45ab09e0e80"
 
 [[package]]
 name = "csaps"
-version = "0.4.0"
+version = "0.5.0"
 dependencies = [
  "almost",
  "approx 0.5.1",
@@ -96,9 +96,9 @@ checksum = "d231dfb89cfffdbc30e7fc41579ed6066ad03abda9e567ccafae602b97ec5024"
 
 [[package]]
 name = "itertools"
-version = "0.13.0"
+version = "0.14.0"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "413ee7dfc52ee1a4949ceeb7dbc8a33f2d6c088194d9f922fb8318faf1f01186"
+checksum = "2b192c782037fadd9cfa75548310488aabdbf3d2da73885b31bd0abd03351285"
 dependencies = [
  "either",
 ]
@@ -127,9 +127,9 @@ dependencies = [
 
 [[package]]
 name = "ndarray"
-version = "0.16.1"
+version = "0.17.2"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "882ed72dce9365842bf196bdeedf5055305f11fc8c03dee7bb0194a6cad34841"
+checksum = "520080814a7a6b4a6e9070823bb24b4531daac8c4627e08ba5de8c5ef2f2752d"
 dependencies = [
  "approx 0.5.1",
  "matrixmultiply",
@@ -256,9 +256,9 @@ checksum = "3c5e1a9a646d36c3599cd173a41282daf47c44583ad367b8e6837255952e5c67"
 
 [[package]]
 name = "sprs"
-version = "0.11.2"
+version = "0.11.4"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "704ef26d974e8a452313ed629828cd9d4e4fa34667ca1ad9d6b1fffa43c6e166"
+checksum = "6dca58a33be2188d4edc71534f8bafa826e787cc28ca1c47f31be3423f0d6e55"
 dependencies = [
  "alga",
  "ndarray",
@@ -281,9 +281,9 @@ dependencies = [
 
 [[package]]
 name = "syn"
-version = "2.0.79"
+version = "2.0.87"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "89132cd0bf050864e1d38dc3bbc07a0eb8e7530af26344d3d2bbbef83499f590"
+checksum = "25aa4ce346d03a6dcd68dd8b4010bcb74e54e62c90c573f394c46eae99aba32d"
 dependencies = [
  "proc-macro2",
  "quote",
@@ -292,18 +292,18 @@ dependencies = [
 
 [[package]]
 name = "thiserror"
-version = "1.0.64"
+version = "2.0.18"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "d50af8abc119fb8bb6dbabcfa89656f46f84aa0ac7688088608076ad2b459a84"
+checksum = "4288b5bcbc7920c07a1149a35cf9590a2aa808e0bc1eafaade0b80947865fbc4"
 dependencies = [
  "thiserror-impl",
 ]
 
 [[package]]
 name = "thiserror-impl"
-version = "1.0.64"
+version = "2.0.18"
 source = "registry+https://github.com/rust-lang/crates.io-index"
-checksum = "08904e7672f5eb876eaaf87e0ce17857500934f4981c4a0ab2b4aa98baac7fc3"
+checksum = "ebc4ee7f67670e9b64d05fa4253e753e016c6c95ff35b89b7941d6b856dec1d5"
 dependencies = [
  "proc-macro2",
  "quote",
diff --git a/Cargo.toml b/Cargo.toml
index d04f617..47babe3 100644
--- a/Cargo.toml
+++ b/Cargo.toml
@@ -1,6 +1,6 @@
 [package]
 name = "csaps"
-version = "0.4.0"
+version = "0.5.0"
 authors = ["Eugene Prilepin <esp.home@gmail.com>"]
 
 description = "Cubic spline approximation (smoothing)"
@@ -17,19 +17,18 @@ edition = "2021"
 
 
 [badges]
-travis-ci = { repository = "espdev/csaps-rs", branch = "master" }
 coveralls = { repository = "espdev/csaps-rs", branch = "master", service = "github" }
 
 
 [dependencies]
 num-traits = "0.2.14"
-ndarray = "0.16.1"
+ndarray = "0.17.2"
 sprs-ldl = "0.10.0"
-sprs = "0.11.2"
+sprs = "0.11.4"
 almost = "0.2.0"
-itertools = "0.13.0"
-thiserror = "1.0.30" 
+itertools = "0.14.0"
+thiserror = "2.0.18"
 
 [dev-dependencies]
 approx = "0.5.1"
-ndarray = {version = "0.16.1", features = ["approx"]}
+ndarray = {version = "0.17.2", features = ["approx"]}
diff --git a/src/ndarrayext.rs b/src/ndarrayext.rs
index a29ad01..f93f6f4 100644
--- a/src/ndarrayext.rs
+++ b/src/ndarrayext.rs
@@ -1,5 +1,5 @@
 use itertools::Itertools;
-use ndarray::{prelude::*, IntoDimension, Slice};
+use ndarray::{prelude::*, IntoDimension, Order, Slice};
 
 use crate::{
     util::dim_from_vec,
@@ -7,6 +7,19 @@ use crate::{
     Real, Result,
 };
 
+pub(crate) fn reshape_order<T, D>(data: &ArrayView<'_, T, D>) -> Order
+where
+    D: Dimension,
+{
+    if data.is_standard_layout() {
+        Order::RowMajor
+    } else if data.ndim() > 1 && data.raw_view().reversed_axes().is_standard_layout() {
+        Order::ColumnMajor
+    } else {
+        Order::RowMajor
+    }
+}
+
 pub fn diff<'a, T: 'a, D, V>(data: V, axis: Option<Axis>) -> Array<T, D>
 where
     T: Real<T>,
@@ -43,7 +56,10 @@ where
     let axis_size = shape[axis.0];
     let new_shape = [numel / axis_size, axis_size];
 
-    match data_view.permuted_axes(axes).into_shape(new_shape) {
+    let data_view = data_view.permuted_axes(axes);
+    let order = reshape_order(&data_view);
+
+    match data_view.into_shape_with_order((new_shape, order)) {
         Ok(view_2d) => Ok(view_2d),
         Err(error) => Err(ReshapeTo2d {
             input_shape: shape,
@@ -62,7 +78,9 @@ where
     let shape = data.shape().to_vec();
     let new_shape = [shape[0..(ndim - 1)].iter().product(), shape[ndim - 1]];
 
-    match data.into_shape(new_shape) {
+    let order = reshape_order(&data);
+
+    match data.into_shape_with_order((new_shape, order)) {
         Ok(data_2d) => Ok(data_2d),
         Err(error) => Err(ReshapeTo2d {
             input_shape: shape,
@@ -93,7 +111,9 @@ where
     let new_shape: D = dim_from_vec(ndim, new_shape_vec.clone());
     let data_view = data.into();
 
-    match data_view.into_shape(new_shape) {
+    let order = reshape_order(&data_view);
+
+    match data_view.into_shape_with_order((new_shape, order)) {
         Ok(view_nd) => {
             let mut axes_tmp: Vec<usize> = (0..ndim).collect();
             let end_axis = axes_tmp.pop().unwrap();
diff --git a/src/ndg/evaluate.rs b/src/ndg/evaluate.rs
index 63f493c..b92d05b 100644
--- a/src/ndg/evaluate.rs
+++ b/src/ndg/evaluate.rs
@@ -1,6 +1,10 @@
 use ndarray::{Array, ArrayView1, Dimension};
 
-use crate::{ndarrayext::to_2d_simple, util::dim_from_vec, NdSpline, Real};
+use crate::{
+    ndarrayext::{reshape_order, to_2d_simple},
+    util::dim_from_vec,
+    NdSpline, Real,
+};
 
 use super::{util::permute_axes, GridCubicSmoothingSpline, NdGridSpline};
 
@@ -36,8 +40,10 @@ where
                 coeffs_shape[ndim_m1] = xi_ax.len();
                 let shape: D = dim_from_vec(self.ndim, coeffs_shape);
 
+                let order = reshape_order(&coeffs_2d.view());
+
                 coeffs_2d
-                    .into_shape(shape)
+                    .into_shape_with_order((shape, order))
                     .unwrap()
                     .permuted_axes(permuted_axes.clone())
                     .to_owned()
diff --git a/src/ndg/make.rs b/src/ndg/make.rs
index 99e8e4d..f0a1edb 100644
--- a/src/ndg/make.rs
+++ b/src/ndg/make.rs
@@ -1,7 +1,10 @@
 use ndarray::Dimension;
 
 use crate::util::dim_from_vec;
-use crate::{ndarrayext::to_2d_simple, CubicSmoothingSpline, Real, RealRef, Result};
+use crate::{
+    ndarrayext::{reshape_order, to_2d_simple},
+    CubicSmoothingSpline, Real, RealRef, Result,
+};
 
 use super::{util::permute_axes, GridCubicSmoothingSpline, NdGridSpline};
 
@@ -54,9 +57,11 @@ where
                 coeffs_shape[ndim_m1] = spline.pieces() * spline.order();
                 let new_shape: D = dim_from_vec(ndim, coeffs_shape);
 
-                spline
-                    .coeffs()
-                    .into_shape(new_shape)
+                let coeffs = spline.coeffs();
+                let order = reshape_order(&coeffs);
+
+                coeffs
+                    .into_shape_with_order((new_shape, order))
                     .unwrap()
                     .permuted_axes(permuted_axes.clone())
                     .to_owned()
diff --git a/src/umv/evaluate.rs b/src/umv/evaluate.rs
index 293b02c..5f92516 100644
--- a/src/umv/evaluate.rs
+++ b/src/umv/evaluate.rs
@@ -44,14 +44,15 @@ where
         // Returns NxM array of coeffs values for given 1xM indices array
         // where N is ndim and M is the size of xi
         let get_indexed_coeffs = |inds: &Array1<usize>| {
-            // Returns Nx1 2-d array of coeffs by given index
-            let coeffs_by_index = |&index| coeffs.slice(s![.., index]).insert_axis(Axis(1));
+            let mut indexed_coeffs = Array2::<T>::zeros((coeffs.nrows(), inds.len()));
 
-            // Get the M-sized vector of coeffs values Nx1 arrays
-            // for all dimensions for 1xM indices array
-            let indexed_coeffs: Vec<_> = inds.iter().map(coeffs_by_index).collect();
+            for (col, &index) in inds.iter().enumerate() {
+                indexed_coeffs
+                    .column_mut(col)
+                    .assign(&coeffs.slice(s![.., index]));
+            }
 
-            concatenate(Axis(1), &indexed_coeffs).unwrap()
+            indexed_coeffs
         };
 
         // Vectorized computing the spline pieces (polynoms) on the given data sites