Solving double pendulum with runge kutta

I have four first order differential equations (see attachement diff.jpg).

where
g= 9.8m/s
l1 = 0.216 m
l2 = 0.241 m
Δθ = θ1-θ2
β = g/l1
α = l2/l1

μ = 1 + m1/m2 = 2.22

// Initial values
t = 0.0;
theta1 = (141.7/180)*pi;
theta2 = (169.8/180)*pi;
p = 0.0;
q = 0.0;
the time step size , h=0.01s

please see the second order runge kutta funtion I created in C++ and advice if the method was correct

void RungeKutta(double *t, double *theta1, double *p, double *theta2, double *q, double h)
{
double tVal, theta1Val, pVal, theta2Val, qVal;
double k1, n1, k2, n2, k, n;
double K1, N1, K2, N2, K, N;

tVal = *t;
theta1Val = *theta1;
pVal = *p;
theta2Val = *theta2;
qVal = *q;

//cout << endl << tVal << " " << theta1Val << " " << pVal << " " << theta2Val << " " << qVal <<endl;

k1 = h*f(tVal, theta1Val, pVal, theta2Val, qVal);
n1 = h*g(tVal, theta1Val, pVal, theta2Val, qVal);
K1 = h*F(tVal, theta1Val, pVal, theta2Val, qVal);
N1 = h*G(tVal, theta1Val, pVal, theta2Val, qVal);


k2 = h*f(tVal+h, theta1Val+k1, pVal+n1, theta2Val+K1, qVal+N1);
n2 = h*g(tVal+h, theta1Val+k1, pVal+n1, theta2Val+K1, qVal+N1);
K2 = h*F(tVal+h, theta1Val+k1, pVal+n1, theta2Val+K1, qVal+N1);
N2 = h*G(tVal+h, theta1Val+k1, pVal+n1, theta2Val+K1, qVal+N1);

//cout << k2 << " " << n2 << " "<< K2 << " "<< N2 << " " <<endl;

k = (k1 + k2)/2;
n = (n1 + n2)/2;
K I have four first order differential equations (see attachement diff.jpg).

where
g= 9.8m/s
l1 = 0.216 m
l2 = 0.241 m

B = g/l1
alpha = l2/l1

u = 1 + m1/m2 = 2.22

// Initial values
t = 0.0;
theta1 = (141.7/180)*pi;
theta2 = (169.8/180)*pi;
p = 0.0;
q = 0.0;
the time step size , h=0.01s

please see the second order runge kutta funtion I created in C++ and advice if the method was correct

void RungeKutta(double *t, double *theta1, double *p, double *theta2, double *q, double h)
{
double tVal, theta1Val, pVal, theta2Val, qVal;
double k1, n1, k2, n2, k, n;
double K1, N1, K2, N2, K, N;

tVal = *t;
theta1Val = *theta1;
pVal = *p;
theta2Val = *theta2;
qVal = *q;

//cout << endl << tVal << " " << theta1Val << " " << pVal << " " << theta2Val << " " << qVal <<endl;

k1 = h*f(tVal, theta1Val, pVal, theta2Val, qVal);
n1 = h*g(tVal, theta1Val, pVal, theta2Val, qVal);
K1 = h*F(tVal, theta1Val, pVal, theta2Val, qVal);
N1 = h*G(tVal, theta1Val, pVal, theta2Val, qVal);


k2 = h*f(tVal+h, theta1Val+k1, pVal+n1, theta2Val+K1, qVal+N1);
n2 = h*g(tVal+h, theta1Val+k1, pVal+n1, theta2Val+K1, qVal+N1);
K2 = h*F(tVal+h, theta1Val+k1, pVal+n1, theta2Val+K1, qVal+N1);
N2 = h*G(tVal+h, theta1Val+k1, pVal+n1, theta2Val+K1, qVal+N1);

//cout << k2 << " " << n2 << " "<< K2 << " "<< N2 << " " <<endl;

k = (k1 + k2)/2;
n = (n1 + n2)/2;
K = (K1 + K2)/2;
N = (N1 + N2)/2;

// cout << endl << k << " " << n << " "<< K << " "<< N << " " <<endl;

tVal += h; theta1Val += k; pVal += n; theta2Val += K; qVal += N;

*t = tVal;
*theta1 = theta1Val;
*p = pVal;
*theta2 = theta2Val;
*q = qVal;

} = (K1 + K2)/2;
N = (N1 + N2)/2;

// cout << endl << k << " " << n << " "<< K << " "<< N << " " <<endl;

tVal += h; theta1Val += k; pVal += n; theta2Val += K; qVal += N;

*t = tVal;
*theta1 = theta1Val;
*p = pVal;
*theta2 = theta2Val;
*q = qVal;

}