Something like this?
clear; close all; clc;
% Generate example data
[X, Y] = meshgrid(-5:0.1:5);
Z = peaks(X, Y);
C = meshgrid(linspace(0.5, 2, size(X,1))); % Fourth array with random values
% Create surface plot
surf(X, Y, Z, C, 'EdgeColor', 'none');
colorbar
clim([min(C(:)), max(C(:))])
CLIM = clim;
n_color = 100;
% finding the corresponding index for specific values like 1 and 1.5
indfor1 = cindex(1, CLIM, n_color);
indfor1p5 = cindex(1.5, CLIM, n_color);
cmap = colormap(parula(n_color));
cmap(1:indfor1,:) = repmat([0,0,1], indfor1, 1);
% transition from [1,1,0] (yellow) to [1, 0, 0] (red)
cYellow = [1, 1, 0]; cRed = [1, 0, 0];
cmap(indfor1:indfor1p5,:) = ...
[linspace(cYellow(1), cRed(1), indfor1p5-indfor1+1)', ...
linspace(cYellow(2), cRed(2), indfor1p5-indfor1+1)', ...
linspace(cYellow(3), cRed(3), indfor1p5-indfor1+1)'];
cmap(indfor1p5+1:end,:) = repmat([1,0,0], n_color-indfor1p5, 1);
colormap(cmap)
function ind = cindex(val, CLIM, n_color)
% To get the index of colormap for a specific value, val
% CLIM: colormap limits
% n_color: the number of colormap
ind = (val - CLIM(1)) / diff(CLIM) * (n_color-1) + 1;
end
