# Using Matlab to solve second order ODE's

Vuldoraq

## Homework Statement

A projectile is ﬁred with an initial speed of 1000 m/s at an angle α from the horizontal. A
good model for the path uses an air-resistance force proportional to the square of the speed.
This leads to the following equations for the acceleration in the x and y directions:

$$\ddot{x}=-c*\dot{x}*(\dot{x}^{2}+\dot{y}^{2})^{1/2}$$
$$\ddot{y}=-c*\dot{y}*(\dot{x}^{2}+\dot{y}^{2})^{1/2}-g$$

where an over dot signiﬁes a derivative with respect to time, and the constant c can be taken to be c = 0.005. The initial conditions are:

$$x(0) =y(0) =0\ and\ (\dot{x(0)},\dot{y(0)}) =1000(cos\alpha,sin\alpha)$$

(b) Convert the above system of equations to 4 ﬁrst-order equations.
(c) Solve the system of equations over the time interval [0,30] (with e.g. the RK4 method)
for each launch angle α =10:10:80 degrees (convert to radians!). Plot the resulting
trajectories.

None

## The Attempt at a Solution

I transformed the second order ode's to four first order;

$$\frac{dx}{dt}=x_{1}$$

$$\frac{dx_{1}}{dt}=-c*x_{1}*(x_{1}^{2}+y_{1}^{2})^{1/2}$$

$$\frac{dy}{dt}=y_{1}$$

$$\frac{dy_{1}}{dt}=-c*y_{1}*(x_{1}^{2}+y_{1}^{2})^{1/2}-g$$

And wrote this program to represent them,

function [dx]=secondorder(t,v)
c=0.005;
g=9.8;

x=v(1);
y=v(2);

dx(1)=x;
dx(2)=-c.*x.*sqrt(x.^2+y.^2);
dx(3)=y;
dx(4)=-c.*y.*sqrt(x.^2+y.^2)-g;

Then fed this program into this fourth order Runge-Kutte program,

function [sol]=rk4_syst(fcn,a,b,y0,N) % fcn now outputs the vector of derivs.

h=(b-a)/N; % and y0 is the vector if initial values
y(1,:) = y0; % copy into the first vector-point
x=a+(0:N)*h;

for k =1:N % loop over steps as before
k1 = feval(fcn,x(k),y(k, ; % for each step, compute
k2 = feval(fcn,x(k)+h/2,y(k,:)+h*(k1)/2); % the slopes k1,...,k4
k3 = feval(fcn,x(k)+h/2,y(k,:)+h*k2/2); % for all the equations i
k4 = feval(fcn,x(k)+h,y(k,:)+h*k3); % by using vector-notation
y(k+1,:)=y(k,:)+h*(k1+2*(k2+k3)+k4)/6; % then take vector step
end % (simultaneously for all i)

sol=[x',y];

But everytime I try to solve my system I get silly results, for instance the y part goes way to high or the x part stays fixed! I've spent ages trying to puzzle it out to no avail. Please could someone give me a hand? I know it's an involved question, but I've run of palces to turn. 