I'd probably concatenate the x position and depth of every salinity measurement, then grid them up with gridfit. Here's an example that you can run for four CTD casts.
% x locations along a transect:
x = [0 10 12 18];
% 30 or so random depths from 0 to 600ish dbar:
d1 = sort(500*rand(30,1));
d2 = sort(540*rand(25,1));
d3 = sort(585*rand(35,1));
d4 = sort(510*rand(28,1));
% Salinity values at each measurement:
s1 = 33 + (.002*d1).^3;
s2 = 33.35 + (.002*d2).^3;
s3 = 33.3 + (.002*d3).^3;
s4 = 33.2 + (.002*d4).^3;
% Plot the salinity measurements:
scatter(x(1)*ones(size(s1)),d1,40,s1,'filled');
hold on
scatter(x(2)*ones(size(s2)),d2,40,s2,'filled');
scatter(x(3)*ones(size(s3)),d3,40,s3,'filled');
scatter(x(4)*ones(size(s4)),d4,40,s4,'filled');
% Make a colorbar:
cb = colorbar;
% Set colormap:
cmocean haline
% Label axes:
xlabel('distance along transect (km)')
ylabel('pressure (dbar)')
ylabel(cb,'salinity')
% reverse order of y axes:
axis ij tight
Now to fill in the space between each measurement the way ODV does we need to interpolate to get a value at each pixel on the screen. There are different ways of interpolating and each method can have drawbacks. With interp1 you can interpolate each profile to common depth values and then use interp1 again to interpolate those interpolated values to distances along the transect, but that will leave some big gaps in the data (and such a method require more manual steps on your part than necessary).
Alternatively, you can use the built in functions griddata or scatteredInterpolant, but both tend to produce nonphysical artifacts like striping or voronoi type patterns. Thus, I prefer the simple, elegant method of using John D'Errico's gridfit like this:
% x values corresponding to each measurement:
xi = [repmat(x(1),size(d1));repmat(x(2),size(d2));repmat(x(3),size(d3));repmat(x(4),size(d4))];
% All depths:
di = [d1;d2;d3;d4];
% All salinity measurments:
si = [s1;s2;s3;s4];
% Define some common distances along transect:
X = 0:0.1:20;
% Define some common depths:
D = 0:5:600;
S = gridfit(xi,di,si,X,D);
p = pcolor(X,D,S) ;
shading interp
uistack(p,'bottom')
% Clarify the location of each measurement:
plot(xi,di,'.','color',0.5*[1 1 1])
% Plot bathymetry:
p = patch([x 20 0],[max(d1) max(d2) max(d3) max(d4) 600 600],'k');
set(p,'facecolor',0.5*[1 1 1])