I hope this will help you:
clc;
clear variables;
close all;
% Parameters of the model
S = rand(1);
M = 10;
B = 2;
g = 9.81;
a = deg2rad(35)*(1+0.2*S);
V0 = 28*(1+0.2*S);
% Dynamic model
% XDot = f(t,X), with X = [Vx,Vy,X,Y]
F = @(t,X) [-B/M*X(1); -g-B/M*X(2); X(1); X(2)];
% Initial conditions
X0 = [V0*cos(a); V0*sin(a); 0; 0];
% Numerical solution
[T,X] = ode45(F,[0,4],X0);
% Plots
subplot(2,2,1); plot(T,X(:,1)); title('Vx vs. Time'); grid on;
ylabel('Vel. [m/s]'); xlabel('Time [s]');
subplot(2,2,2); plot(T,X(:,2)); title('Vy vs. Time'); grid on;
ylabel('Vel. [m/s]'); xlabel('Time [s]');
subplot(2,2,3); plot(T,X(:,3)); title('X vs. Time'); grid on;
ylabel('Pos. [m]'); xlabel('Time [s]');
subplot(2,2,4); plot(T,X(:,4)); title('Y vs. Time'); grid on;
ylabel('Pos. [m]'); xlabel('Time [s]');
