From db207bc9f48278ca4def15bf90591bfa44d0801d Mon Sep 17 00:00:00 2001
From: Behnam Hashemi <hoseynhashemi@gmail.com>
Date: Fri, 15 Sep 2023 00:15:22 +0100
Subject: [PATCH] fix-chebfun3-feval-matrix-case

Addresses issue #2436
---
 @chebfun3/feval.m           |  8 ++++----
 tests/chebfun3/test_feval.m | 18 ++++++++++++++----
 2 files changed, 18 insertions(+), 8 deletions(-)

diff --git a/@chebfun3/feval.m b/@chebfun3/feval.m
index f3a961db79..29bd2ddced 100644
--- a/@chebfun3/feval.m
+++ b/@chebfun3/feval.m
@@ -262,7 +262,7 @@
         max(max(max(abs(bsxfun(@minus, yeval, yeval(:, 1)))))) == 0  && ... 
         max(max(max(abs(bsxfun(@minus, zeval, zeval(1, 1)))))) == 0 )
     
-    xeval = xeval(1, :)';
+    xeval = xeval(1, :).';
     yeval = yeval(:, 1);
     zeval = zeval(1, 1);
     colsVals = f.cols(xeval);
@@ -277,14 +277,14 @@
         max(max(max(abs(bsxfun(@minus, yeval, yeval(1, 1)))))) == 0  && ... 
         max(max(max(abs(bsxfun(@minus, zeval, zeval(:, 1)))))) == 0 )
     
-    xeval = xeval(1, :)'; 
+    xeval = xeval(1, :).'; 
     yeval = yeval(1, 1); 
     zeval = zeval(:, 1); 
     colsVals = f.cols(xeval); 
     rowsVals = f.rows(yeval);
     tubesVals = f.tubes(zeval);
     out = squeeze(chebfun3.txm(chebfun3.txm(chebfun3.txm(f.core, colsVals, 1),...
-        rowsVals, 2), tubesVals, 3));
+        rowsVals, 2), tubesVals, 3)).';
 
 % Data generated by meshgrid & xx is a multiple of a single scalar    
 elseif ( max(max(max(abs(bsxfun(@minus, xeval, xeval(1, 1)))))) == 0  && ...                     
@@ -298,7 +298,7 @@
     rowsVals = f.rows(yeval);
     tubesVals = f.tubes(zeval);
     out = squeeze(chebfun3.txm(chebfun3.txm(chebfun3.txm(f.core,...
-        colsVals, 1), rowsVals, 2), tubesVals, 3));
+        colsVals, 1), rowsVals, 2), tubesVals, 3)).';
     
 else
     % Inputs are not obtained by ndgrid or meshgrid.
diff --git a/tests/chebfun3/test_feval.m b/tests/chebfun3/test_feval.m
index 3877289d3e..dd2561c8d7 100644
--- a/tests/chebfun3/test_feval.m
+++ b/tests/chebfun3/test_feval.m
@@ -111,8 +111,7 @@
 FF = ff(xx,yy,zz);
 pass(18) = norm(F - FF) < 100*tol;
 
-% y and z are matrix inputs (meshgrid) and x is just a matrix of just one 
-% scalar
+% y and z are matrix inputs (meshgrid) and x is a matrix of just one scalar
 ff = @(x,y,z) sin(pi*(x+y+z));
 f = chebfun3(ff);
 zz = linspace(-1, 1, 100)';
@@ -144,8 +143,7 @@
 FF = ff(xx,yy,zz);
 pass(21) = norm(F - FF) < 100*tol;
 
-% y and z are matrix inputs (ndgrid) and x is just a matrix of just one 
-% scalar
+% y and z are matrix inputs (ndgrid) and x is a matrix of just one scalar
 ff = @(x,y,z) sin(pi*(x+y+z));
 f = chebfun3(ff);
 zz = linspace(-1, 1, 100)';
@@ -350,4 +348,16 @@
 pass(47) = all(size (xx) == size(vals));
 pass(48) = norm(ff(xx,yy,zz) - vals) < tol;
 
+%% Test from #2436: Ensure output is not transposed incorrectly.
+% Follows Test 20 above.
+ff = @(x,y,z) x;
+f = chebfun3(ff);
+xx = [-0.4 0.4; 
+      -0.4 0.4];
+yy = [-0.2 -0.2;
+       0.2  0.2];
+FF = ff(xx,0*xx,yy);
+F = f(xx,0*xx,yy); 
+pass(49) = norm(F - FF) < 100*tol;
+
 end
\ No newline at end of file