couplode = @(t,x) [x(2); x(4)^2 + tan(x(3)); x(4);  x(2)+x(4)+cos(x(3))];
[t,x] = ode45(couplode, [0 0.4999*pi], [0;0;0;0]); 
figure(1)
% Lembrando que x = [y, ydot, z, zdot]
plot(t, x)
grid
str = {'$ y $','$ \dot{y} $', '$ z $', '$ \dot{z} $'};
legend(str, 'Interpreter','latex')
% ou
y=x(:,1);
ydot=x(:,2);
z=x(:,3);
zdot=x(:,4);
figure(2)
plot(t,zdot)
legend({'$ \dot{z} $'}, 'Interpreter','latex')