throw solved by ode45

program to solve throw by using ode45

clc
clear
close all
n=500;  %number of iterations
succesful=0;  %number of succesful throws


for i=1:n
speed = 10 + randn*0.5;
angle = 40 + randn*5;

[Diter,Hiter,Titer] = Throw_calc(speed,angle);

D(i)=Diter; %vector of all distances
H(i)=Hiter; %vector of all heights
T(i)=Titer;

if((Diter>=9.9)&(Diter<=10.1))
    succesful=succesful+1;
end

end

function [D,H,T] = Throw_calc(speed,angle)
t=0;
Ts=0.01; %step size
Tmax=15; %max time of iterating

theta=angle*pi/180; %recalculating angle into rad

X0=[0;cos(theta)*speed;0;sin(theta)*speed];

options = odeset('Events',@odeEvents);

[t,X]=ode45(@MyOde,[t:Ts:Tmax],X0,options);


%out of function
D=X(end,1);
H=max(X(:,3));
T=t(end);
%ploting
hold on


if((D>=9.5)&(D<=10.5))
    plot(X(:,1),X(:,3),'g')
else
     plot(X(:,1),X(:,3),'r')
end

end


function [position,isterminal,direction] = odeEvents(t,X)
  position = X(3); 
  isterminal = 1;  
  direction = 0;   
end

function [Dx]=MyOde(t,X)
 g=9.81;
A=[0 1 0 0;
   0 0 0 0;
   0 0 0 1;
   0 0 0 0];
 B=[0;0;0;-g];
Dx=A*X+B;
end

or rungerkhuta

clear all
close all
clc
 
tSpan = 7.2;
m = 2;
h = 10e-2;
v0 = 50;
phi = pi/4;
X0 = [0; cos(phi)*v0; 0 ;sin(phi)*v0];
n = tSpan/h;
t = 0:h:tSpan;
x = X0;
for i = 2:n+1
 
    a1 = odefcn(t(i-1),x(:,i-1));
    a2 = odefcn(t(i-1),x(:,i-1) + h/2*a1);
    a3 = odefcn(t(i-1),x(:,i-1) + h/2*a2);
    a4 = odefcn(t(i-1),x(:,i-1) + h*a3);
 
 
    x(:,i) = x(:,i-1) + h/6*(a1 + 2*a2 + 2*a3 + a4);
end





t_min = 0;
t_max = 10;
x0 = 1; 
h = 1;  
N = (t_max-t_min)/h; 




x_euler(1) = x0;
t=0;
for i = 1:N

    a = odefcn(t(i),x_euler(i));
    x_euler(i+1) = x_euler(i) + h*a;
    t(i+1) = t(i) + h;

end



x_trapezoidal(1) = x0; % IC
t = 0;
for i = 1:n

    a1 = odefcn(t(i),x_trapezoidal(i));
    x_hat = x_trapezoidal(i) + h*a1;
    a2 = odeFcn_LODR1(t(i)+h,x_hat);

    x_trapezoidal(i+1) = x_trapezoidal(i) + h/2*(a1 + a2);
    t(i+1) = t(i) + h;
end




 function [dX] = odefcn(t,x)
 
g = 9.81;
 
A = [0 1 0 0; 0 0 0 0; 0 0 0 1; 0 0 0 0];
B = [0; 0; 0; -g];
dX = A*x + B;
 
end
 
function [dx] = odefcn(t,x)
    dx = -0.5*x; 
end