Here are two related methods. The first way just does a smoothing spline on the new function that you have defined:
c0 = 1800;
g = fit(z,(1./c.^2 - 1/c0^2),'smoothingspline');
I1 = integrate(g,50,30)
Because of the wonky syntax for the integrate function, the limits 30 and 50 have been reversed so that you get a positive result.
Since sound speed is the primary data, you may feel that it's better to fit sound speed with a smoothing spline rather than some derived quantity. In that case
I2 = integral(@(x) (1./f(x).^2 - 1/c0^2)', 30, 50)
works (note that a transpose sign ' is necessary). Since the original function c is so good, the two methods agree to within a couple of parts per million.