If you are okay with using the built-in function, then use ode45, which is an implementation of the RK (4,5) algorithm.
data = load('angular.velocity.data.txt');
v1 = data(:,1); % First column
v2 = data(:,2); % Second column
v3 = data(:,3); % Third column
time = 0:0.016:0.016*(numel(v1)-1);
IC = 0; % initial displacement is 0
[~,x1] = ode45(@(t,x) odeFun(t, x, time, v1), time, 0);
[~,x2] = ode45(@(t,x) odeFun(t, x, time, v2), time, 0);
[~,x3] = ode45(@(t,x) odeFun(t, x, time, v3), time, 0);
%% plotting
subplot(1,3,1)
plot(t,v1,t,x1);
legend({'Velocity', 'displacement'})
subplot(1,3,2)
plot(t,v2,t,x2);
legend({'Velocity', 'displacement'})
subplot(1,3,3)
plot(t,v3,t,x3);
legend({'Velocity', 'displacement'})
function v = odeFun(t, x, time, vx)
% v = interp1(time, vx, t, 'linear'); % linear interpolation
% v = interp1(time, vx, t, 'pchip'); % pchip interpolation
v = interp1(time, vx, t, 'makima'); % makima interpolation
end
