How to vary launching angle until projectile lands at a certain spot?

Hey,

Basically i have to write a program to plot a projectiles trajectory which flies to a certain spot and lands there, it must include drag coefficients, i have already sorted part of it out but i'm stuck on the bit where i need to create something to make it alter the launch angle until it hits the right displacement in the x-axis at the time of landing, any help would be greatly appreciated. Here is my code;

function [introduce_velocity]=trajectory28(Vi)

disp('If there is any error please comfirm that the variable values have been input')         %Reminder to check if the values have been input by the user

%-------------------------------INITIAL VALUES--------------------------------------%

s=10000;            %distance to target
%Vi=600;            %launch speed
m=50;               %projectile mass
tp=15;               %parachute opening time
g=9.81;             %gravitational acceleration
A=0.01;             %projectile area
Cd=0.4;             %projectile drag coefficient
Ap=0.05;            %parachute area
Cdp=1.2;            %parachute drag coefficient
roi=1.207;          %air density at ground level
dt=0.1;             %insignificant increment in time
launch_angle=pi/4;  %launch angle
d_angle=0.1;
x=0;                %initial displacement in the x-axis
y=0;                %initial displacement in the y-axis
t=0;                %time from where calculations start
rho=1.207;          %density at ground level but also because in the formula we need a roi value so this one is the variable with the iteration

%-----------------------EQUIVALENT INITIAL VALUES IN EACH AXIS------------------------%

%while x(I+1)<10000
 %launch_angle=launch_angle+d_angle

vx=Vi*cos(launch_angle);            %initial velocity in the x-axis because Vi is a value the user inputs
vy=Vi*sin(launch_angle);            %initial velocity in the y-axis because Vi is a value the user inputs

vv=sqrt(vx^2+vy^2);                 %modulus of initial velocity changed name for code iteration purposes

D=(roi*(vv^2)*Cd*A)/(-2);           %initial drag value

theta=pi/4;                         %angle at which the iteration starts, obviously is the same as the launch_angle

hold on                             %to keep all the values throughout the iteration for plotting purposes

   %---------------------ITERATION FOR ALL THE VALUES FROM TIME '0' TO PARACHUTE OPENING---------------------%

for I=1:10*tp

      D(I+1)=(rho(I)*(vv(I)^2)*Cd*A)/(-2);        %Drag at every point until parachute deployment

      ax(I+1)=(D(I)*cos(theta(I)))/m;             %Acceleration in x-axis at any point until parachute deployment
      ay(I+1)=(D(I)*sin(theta(I))-m*g)/m;         %Acceleration in y-axis at any point until parachute deployment

      aa(I+1)=sqrt((ax(I+1))^2+(ay(I+1))^2);      %Modulus of acceleration at any point until parachute deployment

      vx(I+1)=vx(I)+dt*ax(I);                     %Velocity in the x-axis at any point until parachute deployment
      vy(I+1)=vy(I)+dt*ay(I);                     %Velocity in the y-axis at any point until parachute deployment

      x(I+1)=x(I)+dt*vx(I);                       %Displacement in the x-axis at any point until parachute deployment
      y(I+1)=y(I)+dt*vy(I);                       %Displacement in the y-axis at any point until parachute deployment

      max_displacementx_of_vector=max(x(I+1));    %Syntax to calculate maximum value of displacement vector, for plotting purposes
      max_displacementy_of_vector=max(y(I+1));    %Syntax to calculate maximum value of displacement vector, for plotting purposes

      vv(I+1)=sqrt(vx(I)^2+vy(I)^2);              %Modulus of velocity at any point until parachute deployment

      rho(I+1)=roi*(1-2.333*10^(-5)*y(I+1))^5;    %Changing density as a result of height increase

      if y(I+1)>42850                             %If statement to break the loop in case of a height increase of 42850 due to that formula no longer being valid
          break
      end

      theta(I+1)=atan(vy(I)/vx(I));               %By using the inverse tangent of velocity in the y-axis divided by the x-axis we obtain the new angle for the drag formula
      degrees=(theta(I)*360)/2*pi;                %Formula to convert the theta angle from radians to degrees for ease of reading purposes

      t=t+dt;                                     %By using this formula we obtain the desired iterations

      %----------------------------PLOTTING OF THE TRAJECTORY PRIOR TO PARACHUTE------------------------------%
      %plot(a,'g*',b)
      plot(t,max_displacementx_of_vector,'g*')
      plot(t,max_displacementy_of_vector,'b+')
      %plot(t,vx,'g*')
      %plot(t,vy,'yo')
      %plot(t,ax,'p*')
      %plot(t,ay,'w*')
  end

for I2=10*tp:1500

      D(I2+1)=(rho(I2)*(vv(I2)^2)*Cdp*Ap)/(-2);

      ax(I2+1)=(D(I2)*cos(theta(I2)))/m; %acceleration in x-axis
      ay(I2+1)=(D(I2)*sin(theta(I2))-m*g)/m; %acceleration in y-axis

      aa(I2+1)=sqrt((ax(I2+1))^2+(ay(I2+1))^2);

      vx(I2+1)=vx(I2)+dt*ax(I2);%calculate new velocity to be able to repeat step to column above;
      vy(I2+1)=vy(I2)+dt*ay(I2);

      x(I2+1)=x(I2)+dt*vx(I2);%changing displacement based onn initial values
      y(I2+1)=y(I2)+dt*vy(I2);

      a=max(x(I2+1));
      b=max(y(I2+1));

      vv(I2+1)=sqrt(vx(I2)^2+vy(I2)^2);

      rho(I2+1)=roi*(1-2.333*10^(-5)*y(I2+1))^5;

      if y(I2+1)>42850
          break
      end

      theta(I2+1)=atan(vy(I2)/vx(I2));
      degrees=(theta(I2)*360)/2*pi;
      t=t+dt;
      if b<=0
          break
      end
      %plot(a,'y*',b,'r+')
      plot(t,a,'y*')
      plot(t,b,'r+')
      %plot(t,vx,'g*')
      %plot(t,vy,'yo')
      %plot(t,ax,'p*')
      %plot(t,ay,'w*')
  end
  max_iteration=fix(max(I2));
  max_y=max(b)
  %end

speed_at_impact=vv(max_iteration)
acceleration_at_impact=aa(max_iteration)

end