Runge-Kutta algorithm for a given ODE system consider the system given by:
$$x'_{1}=9x_{1}+24x_{2}+5\cos t-\dfrac{1}{3}\sin t$$
$$x'_{2}=-24x_{1}-51x_{2}-9\cos t+\dfrac{1}{3}\sin t$$
with initial values $$x_{1}(0)=\dfrac{4}{3}$$ and $$x_2(0)=\dfrac{2}{3}$$
Implement the RK4 algorithm in MATLAB or JAVA to solve the system from t=0 to t=20. Let h=1/1000,1/25, 1/20,1/15,4/55, compute $x_2(20)$ in each case.
Thanks, i don't even try because i really don´t like numerical analysis nor applied math, so i don´t know how to do that in a programming language, i would sincerely appreciate your help on this, i don´t know how to proceed on that so i would appreciate even a short explanation of what i need to do or how to do that, if you consider this question do not belongs to this stack exchange please tell me and i would retry it from here and put in on a programming or CS forum.
THANKS A LOT
 A: Consider the vector $X=\binom {x_1}{x_2}$ such that $X'=\binom {x'_1}{x_2'}=\binom {9x_{1}+24x_{2}+5\cos t-\dfrac{1}{3}\sin t}{-24x_{1}-51x_{2}-9\cos t+\cfrac{1}{3}\sin t}=\phi(t, X)$.
Now consider two sequences $a_n\sim x_1$ and $b_n\sim x_2$ which implies that $X_n=\binom {a_n}{b_n}$
The $RK4$ method is given by $$\left \{\begin{array}{ll}K_1=\phi(t_n, X_n)\\ K_2=\phi(t_n+\cfrac h2,X_n+\cfrac h2 K_1)\\ K_3=\phi(t_n+\cfrac h2,X_n+\cfrac h2 K_2)\\K_4=\phi(t_n+h,X_n+hk_4)\\ X_{n+1}=X_n+\cfrac h6(K_1+2K_2+2K_3+K_4)\\ X_0=\binom {x_1(0)}{x_2(0)}\end{array}\right.$$
If you can turn this into a code on $MATLAB$ then you're good.
A: To help you get started, here's a Matlab program that solves this system using Euler's method.  You'll need to change the line where $f$ is defined.  (This code is written to be clear, not to be as efficient as possible.)
%This code uses Euler's method to solve the system
% x'(t) = f(x(t),t).
% x(t) is a 2 x 1 column vector whose first component is
% x1(t) and whose second component is x2(t).
% The initial condition is x(0) = x0.

f = @(x,t) [0;0]; % DEFINE f CORRECTLY HERE.

x0 = [4/3 ; 2/3];

h = 1/1000;

t = 0:h:20;
numtVals = length(t);

x = zeros(2,numtVals);

x(:,1) = x0;

for i = 1:(numtVals - 1)

    ti = t(i);
    xi = x(:,i);

    xip1 = xi + h*f(xi,ti);
    x(:,i+1) = xip1;

end

figure('Name','x1 and x2')
plot(x(1,:))
hold on
plot(x(2,:),'color','green')
hold off

