hi there
i am trying to solve a simple problem using matlab. first, a bit of background for you:
a projectile is launched from a canon. after 15 seconds a parachute opens, the aim is to land the projectile exactly 10km away. various constants for drag, mass, parachute area etc are given in the problem. the aim of the excerise is to be able to 'put in' a desired distance at which to land the projectile and have matlab calculate the required angle to fire the canon at

I have written a code so far that uses the runge kutta method to solve the ODE, such that putting in the initial values and an angle will tell you where the projectile will land.

the code so far consists of 3 functions. it is run using the following command in matlab: Solve_Trajectory(0.1,0,0,900,45) where the 45 is the angle of launch and can be varied. the code looks like this:


function 1:

Calculate the acceleration

function [xdd, ydd] = find_acceleration(t, x, y, xd, yd, theta)



P = 1.207; % Air Density at sea level (kg/m^3)

if t <= 15, Apr = 0.25;
else Apr = 1.1; % Projectile frontal area (m^2) (changes at t=15 when parachute opens)
end
   
M = 600; % Mass (kg)
   
g = 9.81; % Gravity

if t <=15, Cd = 0.38;
else Cd= 0.9; % Projectile drag coefficient (changes at t=15 when parachute opens)
end

Velocity = ((xd^2)+(yd^2)); % Veloctiy squared (for drag)

Drag = 0.5*P*Velocity*Apr*Cd; % Drag


% Calculations used to calculate xdd and ydd

xdd = (-Drag*cos(theta))/M;
ydd = ((-Drag*sin(theta))-M*g)/M;


function 2:
% Move the problem one step forward in time using Runge Kutta Method

function [x, xd, xdd, y, yd, ydd, theta] = Runge_kutta(t, x, xd, y, yd, theta, tInc)

%Calculate current acceleration

[xdd, ydd] = find_acceleration(t, x, y, xd, yd, theta);


%Define A,B,C,D, Beta and Delta for both x and y

%Ay and Ax

Ax=(tInc/2)*xdd;
Ay=(tInc/2)*ydd;

%Betay and Betax

Betax=(tInc/2)*(xd+(Ax/2));
Betay=(tInc/2)*(yd+(Ay/2));

[xdd, ydd] = find_acceleration(t+(tInc/2), x+Betax, y+Betay, xd+Ax, yd+Ay, theta);

%By and Bx

Bx=(tInc/2)*xdd;
By=(tInc/2)*ydd;

[xdd, ydd] = find_acceleration(t+(tInc/2), x+Betax, y+Betay, xd+Bx, yd+By, theta);

%Cy and Cx

Cx=(tInc/2)*xdd;
Cy=(tInc/2)*ydd;

%Deltay and Deltax

Deltax=tInc*(xd+Cx);
Deltay=tInc*(yd+Cy);

[xdd, ydd] = find_acceleration(t+(tInc), x+Deltax, y+Deltay, xd+2*Cx, yd+2*Cy, theta);

%Dy and Dx

Dx=(tInc/2)*xdd;
Dy=(tInc/2)*ydd;


%Calculate next set of values

x = x + tInc*(xd+(1/3)*(Ax+Bx+Cx));
xd = xd + (1/3)*(Ax+2*Bx+2*Cx+Dx);

y = y + tInc*(yd+(1/3)*(Ay+By+Cy));
yd = yd + (1/3)*(Ay+2*By+2*Cy+Dy);

theta = atan(yd/xd);

[xdd, ydd] = find_acceleration(t, x, y, xd, yd, theta);


function 3:

% Solve the initial-value ODE and plot the trajectory of the projectile

% inputs
% tInc = time increment (s)
% xAlpha = start horizontal displacement (m)
% yAlpha = start vertical displacement (m)
% Exit_V = exit velocity of the projectile at launch (m/s)
% angle = Angle of launch of the projectile (degrees)

function [t, x, y] = Solve_Trajectory(tInc, xAlpha, yAlpha, Exit_V, angle)

% Set initial conditions

theta(1) = angle*(pi/180); % Convert launch angle into radians
t(1)=0;
x(1) = xAlpha;
xd(1) = Exit_V*cos(theta);
y(1) = yAlpha;
yd(1) = Exit_V*sin(theta);

% Loop through until projectile touches down
i=2;
while y >= 0
    [x(i), xd(i), xdd(i-1), y(i), yd(i), ydd(i-1), theta(i)] = Runge_kutta(t(i-1), x(i-1), xd(i-1), y(i-1), yd(i-1), theta(i-1), tInc);
    t(i)=t(i-1)+tInc;
    i=(i+1);
end

% Plot the results

plot(x, y);
title('Trajectory')
xlabel('Horizontal Distance [m]')
ylabel('Vertical Height [m]')

a = i-1;
x = x(a)


when running Solve_Trajectory(0.1,0,0,900,45) it runs through all of that and plots the trajectory of the projectile, and gives the distance travelled, x.

i appreciate thats a lot to look through, and if you're still reading then i thank you!

i now wish to write a shooting method solver to deduce what angle must be put in to the Solve_Trajectory function to output x=10000

can anyone help me with how to do this? my matlab skill are not great and it's taken me weeks to get this far!

many thanks if anyone can help
:)