Here is your existing code:
%B-spline interpolation
eff96_spline(:,1) = spline(1:numel(eff96(:,1)), eff96(:,1), 1:0.1:numel(eff96(:,1)))';
eff96_spline(:,2) = spline(1:numel(eff96(:,1)), eff96(:,2), 1:0.1:numel(eff96(:,1)))';
eff95_spline(:,1) = spline(1:numel(eff95(:,1)), eff95(:,1), 1:0.1:numel(eff95(:,1)))';
eff95_spline(:,2) = spline(1:numel(eff95(:,1)), eff95(:,2), 1:0.1:numel(eff95(:,1)))';
eff94_spline(:,1) = spline(1:numel(eff94(:,1)), eff94(:,1), 1:0.1:numel(eff94(:,1)))';
eff94_spline(:,2) = spline(1:numel(eff94(:,1)), eff94(:,2), 1:0.1:numel(eff94(:,1)))';
eff90_spline(:,1) = spline(1:numel(eff90(:,1)), eff90(:,1), 1:0.1:numel(eff90(:,1)))';
eff90_spline(:,2) = spline(1:numel(eff90(:,1)), eff90(:,2), 1:0.1:numel(eff90(:,1)))';
eff86_spline(:,1) = spline(1:numel(eff86(:,1)), eff86(:,1), 1:0.1:numel(eff86(:,1)))';
eff86_spline(:,2) = spline(1:numel(eff86(:,1)), eff86(:,2), 1:0.1:numel(eff86(:,1)))';
% creating speed torgue and efficiency variables
lng = length(eff96_spline);
speed(:,1) = eff86_spline(:,1);
speed(:,2) = eff90_spline(:,1);
speed(:,3) = eff94_spline(:,1);
speed(:,4) = eff95_spline(:,1);
speed(:,5) = eff96_spline(:,1);
torgue(:,1) = eff86_spline(:,2);
torgue(:,2) = eff90_spline(:,2);
torgue(:,3) = eff94_spline(:,2);
torgue(:,4) = eff95_spline(:,2);
torgue(:,5) = eff96_spline(:,2);
efficiency = ones(lng,5);
efficiency(:,1) = 86 * ones(1, lng)';
efficiency(:,2) = 90 * ones(1, lng)';
efficiency(:,3) = 94 * ones(1, lng)';
efficiency(:,4) = 95 * ones(1, lng)';
efficiency(:,5) = 96 * ones(1, lng)';
% plot(eff96_spline(:,1),eff96_spline(:,2),'-');
% hold on;
% plot(eff95_spline(:,1),eff95_spline(:,2),'-');
% hold on;
% plot(eff94_spline(:,1),eff94_spline(:,2),'-');
% hold on;
% plot(eff90_spline(:,1),eff90_spline(:,2),'-');
% hold on;
% plot(eff86_spline(:,1),eff86_spline(:,2),'-');
% hold on;
% plot(continoustorgue(:,1),continoustorgue(:,2),'k--')
% hold on;
% plot(peaktorgue(:,1),peaktorgue(:,2), 'k--')
% xlabel('Motor speed [rpm]')
% ylabel('Torgue [Nm]')
% xlim([0,5510])
% ylim([0,250])
% legend('96%','95%','94%','90%','86%', 'continous torgue', 'peak torgue')
grid on;
figure()
% surfc([eff86(:,1) eff90(:,1) eff94(:,1), eff95(:,1),eff96(:,1)], [eff86(:,2) eff90(:,2) eff94(:,2),eff95(:,2),eff96(:,2)], (ones(numel(eff86(:,1)),5).*[86 90 94 95 96]))
surfc([eff86_spline(:,1),eff90_spline(:,1), eff94_spline(:,1), eff95_spline(:,1), eff96_spline(:,1)],[eff86_spline(:,2),eff90_spline(:,2), eff94_spline(:,2), eff95_spline(:,2), eff96_spline(:,2)], (ones(lng,5).*[86 90 94 95 96]), ...
"EdgeColor","none", "FaceAlpha",0.5)
xlabel('Motor Speed (rpm]')
ylabel('Torque [Nm]')
zlabel('Efficiency [%]')
colormap(parula(10))
colorbar
Here is what I would add
hold on;
n = length(eff86_spline); % find the length of the data we're working with
foo = [eff86_spline; eff90_spline; eff94_spline; eff95_spline; eff96_spline]; % combine all the speed and torque data into single array
bar = repelem([86; 90; 94; 95; 96], n); % create an array for the efficiency values
blah = [foo, bar]; % combine the speed, torque, and efficiency
blah = unique(blah, "rows"); % remove duplicate data
F = scatteredInterpolant(blah(:,1), blah(:,2), blah(:,3)); % create a scatteredInterpolant object from your data
xq = 2500; % pick a speed
yq = 50; % pick a torque
zr = F(xq, yq); % find the resulting interpolated efficiency
scatter3(xq, yq, zr, "red", "filled") % plot the resulting point
Here is the resulting graph (note I changed the transparency of your surface because the point was slightly hidden below the surface):
