alpha = [10 15 20]; Ma = [3.5 5 8 10 15 20 23];
alpha1 = 10:0.1:20; Ma1 = 3.5:0.1:23;
[Ma1, alpha1] = meshgrid(Ma1, alpha1);
CL = readmatrix('simulation.xlsx', 'Sheet', 'Sheet1', 'Range', 'B2:H4');
CL1 = interp2(Ma, alpha, CL, Ma1, alpha1, 'spline');
CD = readmatrix('simulation.xlsx', 'Sheet', 'Sheet1', 'Range', 'B7:H9');
CD1 = interp2(Ma, alpha, CD, Ma1, alpha1, 'spline');
Sm = 0.484;
m = 907.18;
g = 9.80665;
y0_1 = [6500; 0; 0; 0; 0; 80000];
tspan1 = 0:0.02:300;
options = odeset('Events', @odeEventFun);
[t1, y1] = ode45(@(t,y) Fun1(t, y, m, Sm, CL1, CD1, g), tspan1, y0_1, options);
y0_2 = y1(end, :);
tspan2 = t1(end):0.02:2000;
[t2, y2] = ode45(@(t,y) Fun2(t, y, m, Sm, CL1, CD1, g), tspan2, y0_2);
y0_3 = y2(end, :);
tspan3 = 2000:0.02:2300;
options1 = odeset('Events', @odeEventFun1);
[t3, y3] = ode45(@(t,y) Fun3(t, y, m, Sm, CL1, CD1, g), tspan3, y0_3, options1);
t = [t1; t2; t3];
y = [y1; y2; y3];
Vx = y(:,1); Vy = y(:,2); Vz = y(:,3);
x = y(:,4); y_pos = y(:,5); z = y(:,6);
V = sqrt(Vx.^2 + Vy.^2 + Vz.^2);
figure(1)
plot(t, z