The function to be integrated computes these values. You have all required inputs as output of ODE23S already. Use them as inputs to your function. Convert the function such, that it accepts matrices as inputs.
tspan = [0, 10];
y0 = [0, 0];
[T, Y] = ode23s(@fcn, tspan, y0);
Z = fcn(T.', Y.').';
plot(T, Y); hold('on');
plot(T, Z(:, 2), 'co');
function dy = fcn(t, y)
dy = [y(2, :); ...
sin(t)]; % [EDITED, T -> t]
end
