Confused by Runge-Kutta method

Question

I need to solve flight of the shell problem (get XY trajectory using Runge-Kutta method) and at first check if method works on some simple test system.

So here's Runge-Kutta stuff

$k_1 = f(t_n, y_n)$

$k_2 = f(t_n + h/2, y_n + hk_1/2) $

$k_3 = f(t_n+h, y_n - hk_1 + 2hk_2)$

$y_{n+1} = y_n + h(k_1 + 4k_2 + k_3)/6$

where $h$ is step

Here's my test system

$y'_1 = -5y_1 - 10y_2 + 14e^{-x}$

$y'_2 = -10y_1 - 5y_2 + 14e^{-x}$

approximately it's $y1 = y2 = e^{-x}$

So I solve it on $[0;4]$ and checked it (by comparison with $y1 = y2 = e^{-x}$) (btw approximation is very good even with h=0.1). Here's the matlab code, if you interested in

Rnge-Kutta method

  function [ res_y ] = RungeKutta(dim, size, grid, step, f1,f2,y1, y2)

    k1=zeros(dim);
    k2=zeros(dim);
    k3=zeros(dim);

    h = step;

    res_y(1,1) = y1;
    res_y(2,1) = y2;

    for i=1: size

       k1(1)= f1(grid(i),y1,y2);
       k1(2)= f2(grid(i),y1,y2);

       k2(1)= f1(grid(i)+h/2, y1+h*k1(1)/2, y2+h*k1(2)/2);
       k2(2)= f2(grid(i)+h/2, y1+h*k1(1)/2, y2+h*k1(2)/2);

       k3(1)= f1(grid(i)+h, y1-h*k1(1)+2*h*k2(1), y2-h*k1(2)+2*h*k2(2));
       k3(2)= f2(grid(i)+h, y1-h*k1(1)+2*h*k2(1), y2-h*k1(2)+2*h*k2(2));

       res_y(1,i+1) = y1 + h*(k1(1) + 4*k2(1) + k3(1))/6;
       res_y(2,i+1) = y1 + h*(k1(2) + 4*k2(2) + k3(2))/6;

       y1 = res_y(1,i+1); 
       y2 = res_y(2,i+1);
    end

 end

Main method

a = 0; b = 4;
h = 0.1; % step
t = a:h:b; %grid
n = 2; 
m = size(t,2);


 hold on;
     plot(t, exp(-t),'b-')
     plot(t, exp(-t),'r--')
 hold off;

y1=1; y2 = 1;

f1_ptr = @f1;% out = -5 * y1 - 10 * y2 + (14)*exp(-x);
f2_ptr = @f2;% out = -10 * y1 - 5 * y2 + (14)*exp(-x);

res_y = RungeKutta(n,m-1,t,h,f1_ptr, f2_ptr,1,1);

hold on;
plot(t,res_y);

hold off;

f1 function

function [ out ] = f1( x, y1, y2, alpha, beta )
if nargin == 3
    alpha = 5;
    beta = 10;
end

out = -alpha * y1 - beta * y2 + (alpha + beta - 1)*exp(-x);

end

f2 function

function [ out ] = f2( x, y1, y2, alpha, beta )
if nargin == 3
    alpha = 5;
    beta = 10;
end
out = -beta * y1 - alpha * y2 + (alpha + beta - 1)*exp(-x);

end

But my main system, for the shell trajectory look like that

$x'' = - \frac{1}{2m} C \rho S \cos(\theta) v^2$

$y'' = - \frac{1}{2m} C \rho S \sin(\theta) v^2 - g$

Where's $v = \sqrt{x'^2 + y'^2}$ and angle $\theta = arctg(y'/x')$ and $C, \rho, S, m, g$ some constants.

Also I've got some initial conditions

$x(0) = 0, y(0) = 0, \theta(0) = 0.6, v(0) = 50 m/s$

So I've stuck here. I don't understand clearly how to solve this system using Runge-Kutta method. I mean, I don't understand how to code it in such situation.

I'm not good at math, I am a coder. Please, colud you help me?

UPD Thanks to Tony I think I began to understand problem. But method I wrote doesn't work :c

Here it is

function [ shell_result ] = RungeKuttaShell(step, y1,y2,y3,y4)

    h = step; norma = 1;
    i = 1;
    shell_result(1,1) = y1;
    shell_result(2,1) = y2;
    shell_result(3,1) = y3;
    shell_result(4,1) = y4;

    while (norma > 0) %how much iterations?

       k1(1)= fz(y3);
       k1(2)= fw(y4);
       k1(3)=fx(k1(1), k1(2));
       k1(4)=fy(k1(1), k1(2));

       k2(1)= fz( y3+h*k1(1)/2);
       k2(2)= fw( y4+h*k1(2)/2);
       k2(3)=fx(k2(1),k2(2));
       k2(4)=fy(k2(1),k2(2));

       k3(1)= fz(y3-h*k1(1)+2*h*k2(1));
       k3(2)= fw(y4-h*k1(2)+2*h*k2(2));
       k3(3)=fx(k3(1),k3(2));
       k3(4)=fy(k3(1),k3(2));

       shell_result(1,i+1) = y1 + h*(k1(1) + 4*k2(1) + k3(1))/6;
       shell_result(2,i+1) = y2 + h*(k1(2) + 4*k2(2) + k3(2))/6;
       shell_result(3,i+1) = y3 + h*(k1(3) + 4*k2(3) + k3(3))/6;
       shell_result(4,i+1) = y4 + h*(k1(4) + 4*k2(4) + k3(4))/6;

       y1 = shell_result(1,i+1); 
       y2 = shell_result(2,i+1);
       y3 = shell_result(3,i+1); 
       y4 = shell_result(4,i+1);
       i = i+1;
       norma = y4; %height?
    end

end

functions look like this

function [ output ] = fx( z,w)
    m=15; C=0.2; rho=1.29; S=0.25; g = 9.81;

        output = -((1/2*m)*C*rho*S)*z*sqrt(z*z + w*w);

end

function [ output ] = fy( z,w )
    m=15; C=0.2; rho=1.29; S=0.25; g = 9.81;

    output = -((1/2*m)*C*rho*S)*w*sqrt(z*z + w*w) - g;

end

function [ output_args ] = fw( input_args )

    output_args = input_args;

end

function [ output_args ] = fz( input_args )

     output_args = input_args;

end

And I use it just like that

shell = RungeKuttaShell(0.1,0,0,v0*cos(th0), v0*sin(th0));

Answer 1

The system is $$\begin {cases} x'=z \\ \\y'=w \\ \\ z'=-k\,z\sqrt{z^2+w^2} \\ \\w'=-k\,w\sqrt{z^2+w^2}-g \end {cases}$$with$$k=\frac 1{2m}C\rho S$$The initial values are $x(0)=0\;,\; y(0)=0\;,\;z(0)=v(0)\cos\theta(0)\;,\;w(0)=v(0)\sin\theta(0)\,$.

You must be clever in recoding your procedure.
The system now is $$\begin {cases} y'_1=f_1(t,y_1,y_2,y_3,y_4) \\ \\y'_2=f_2(t,y_1,y_2,y_3,y_4) \\ \\ y'_3=f_3(t,y_1,y_2,y_3,y_4) \\ \\y'_4=f_4(t,y_1,y_2,y_3,y_4) \end {cases}$$The previous one was$$\begin {cases} y'_1=f_1(t,y_1,y_2) \\ \\y'_2=f_2(t,y_1,y_2) \end {cases}$$You must recode completely without any reference to your particular system (this only conditions the choice of the functions).

Now $y_1=x\;$,$\;y_2=y\;$,$\;y_3=z\;$,$\;y_4=w$.

Edit

The choice of step value is not so important if one is not interested in saving resources (of time and space of memory). Do the changes in my comment and try.

If you don't get a range of 176 m, try this

final code

k2(1)= fz( y3+h*k1(3)/2);
k2(2)= fw( y4+h*k1(4)/2);
k2(3)=fx(k2(1),k2(2));
k2(4)=fy(k2(1),k2(2));

k3(1)= fz(y3-h*k1(3)+2*h*k2(3));
k3(2)= fw(y4-h*k1(4)+2*h*k2(4));
k3(3)=fx(k3(1),k3(2));
k3(4)=fy(k3(1),k3(2));