# Runge kutta method Matlab 5th order

How can we moditify the code of Runge Kutta method 2th order from site http://www.computeraideddesignguide.com/euler-vs-runge-kutta-circle-matlab-code/ into the code of Runge Kutta method 5th order?

i tried, but returns back like a spiral. i used Dormand Prince tableau. Can you help to find my mistakes in my code?

 clc;
clear all;
h=0.01;
N=5000;
x=zeros(1,N);
y=zeros(1,N);
x(1)=1; %initial value x(0)
y(1)=0; %initial value y(0)
%iteration with N step
% Forward Euler method
c2=0.2; c3 = 0.3 ;c4=0.8 ;c5=8./9.; c6=1.; c7=1.;
a21 = 0.2; a31 = 3./40.0; a32 = 9.0/40.0; a41 = 44.0/45.0;  a42=56.0/15.0; a43=32.0/9.0;
a51=19372.0/6561.0; a52=25360.0/2187.0; a53=64448.0/6561.0; a54=212.0/729.0;
a61=9017.0/3168.; a62=355.0/33.0; a63=46732.0/5247.0; a64=49.0/176.0; a65=5103.0/18656.0;
a71=35.0/384.0; a73=500.0/1113.0; a74=125.0/192.0 ;a75=2187.0/6784.0; a76=11.0/84.0;
for i=1:N
xk1=y(i);
yk1=-x(i);
xk2 = y(i) + c2*h*(a21*xk1);
yk2 = -(x(i) + c2*h*(a21*yk1));
xk3 = y(i) + c3*h*(a31*xk1 + a32*xk2);
yk3 = -(x(i) + c3*h*(a31*yk1 + a32*yk2));
xk4 = y(i) + c4*h*(a41*xk1 - a42*xk2 + a43*xk3);
yk4 = -(x(i) + c4*h*(a41*yk1 - a42*yk2 + a43*yk3));
xk5 = y(i) + c5*h*(a51*xk1 - a52*xk2 + a53*xk3 - a54*xk4);
yk5 = -(x(i) + c5*h*(a51*yk1 - a52*yk2 + a53*yk3 - a54*yk4 ));
xk6 = y(i) + c6*h*(a61*xk1 - a62*xk2 + a63*xk3 + a64*xk4 - a65*xk5);
yk6 = -(x(i) + c6*h*(a61*yk1 - a62*yk2 + a63*yk3 + a64*yk4 - a65*yk5));
xk7 = y(i) + c7*h*(a71*xk1 + a73*xk3 + a74*xk4 - a75*xk5 + a76*xk6);
yk7 = -(x(i) + c7*h*(a71*yk1 + a73*yk3 + a74*yk4 - a75*yk5 + a76*yk6));
x(i+1) = x(i) + h *((5179.0/57600.)*xk1 + (7571.0/16695.0)*xk3 + (393.0/640.0)*xk4 - (92097.0/339200.0)*xk5 +(187.0/2100.0)*xk6 + (1.0/40)*xk7);
y(i+1) = y(i) + h * ((5179.0/57600.)*yk1 + (7571.0/16695.0)*yk3 + (393.0/640.0)*yk4 - (92097.0/339200.0)*yk5 +(187.0/2100.0)*yk6 + (1.0/40)*yk7);
end
angle=linspace(0,2*pi); % used to plot the unit circle
plot(x,y,'r',sin(angle),cos(angle),'b')
legend('Range Kutta', 'Unit circle')


You need to use the xk1,xk2,... as updates for the x and the yk1,yk2,... as updates for the y. At the moment you have that mixed up, which reduces the order of the method to 1 and you get spirals like with the Euler method. Also, you are using the c values wrong. As the ODE does not depend on time, they should not even be present.
 xk2 = y(i) + h*(a21*yk1);