From 1d847748044ba1dae36b617411a945d3842d1d18 Mon Sep 17 00:00:00 2001
From: danfortunato <dan.fortunato@gmail.com>
Date: Fri, 21 Jul 2023 16:24:48 -0400
Subject: [PATCH] Vectorize chebfun3/feval in the case of unstructured
 evaluation points

---
 @chebfun3/feval.m | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/@chebfun3/feval.m b/@chebfun3/feval.m
index 2931740752..f3a961db79 100644
--- a/@chebfun3/feval.m
+++ b/@chebfun3/feval.m
@@ -113,14 +113,14 @@
     zVec = zVec.';
 end
 
-outVec = zeros(size(xVec));
 colsVals = f.cols(xVec); % So, colsVals(:,j) = feval(f.cols(:,j),xVec).
 rowsVals = f.rows(yVec);
 tubesVals = f.tubes(zVec);
-for i = 1:size(xVec, 1) % TODO: How to vectorize this?
-    outVec(i) = chebfun3.txm(chebfun3.txm(chebfun3.txm(f.core,...
-        colsVals(i,:), 1), rowsVals(i,:), 2), tubesVals(i,:), 3);
-end
+
+txm1 = chebfun3.txm(f.core, colsVals, 1);
+txm2 = pagemtimes(permute(txm1, [3 2 1]), permute(rowsVals, [2 3 1]));
+txm3 = sum(txm2 .* permute(tubesVals, [2 3 1]), 1);
+outVec = reshape(txm3, size(xVec));
 
 if ( (m == 1) && (n > 1) )
     outVec = outVec.';