function [x1,q1,r1] = rk4 (x0,q0,r0,tau,ffun,param);

     [X1,Q1,R1] = deriv(x0,q0,r0,ffun,param);
     [X2,Q2,R2] = deriv(x0+tau/2*X1,q0+tau/2*Q1,r0+tau/2*R1,ffun,param);
     [X3,Q3,R3] = deriv(x0+tau/2*X2,q0+tau/2*Q2,r0+tau/2*R2,ffun,param);
     [X4,Q4,R4] = deriv(x0+tau*X3,q0+tau*Q3,r0+tau*R3,ffun,param);
     
     x1 = x0 + tau/6*(X1 + 2*X2 + 2*X3 + X4);
     q1 = q0 + tau/6*(Q1 + 2*Q2 + 2*Q3 + Q4);
     r1 = r0 + tau/6*(R1 + 2*R2 + 2*R3 + R4);
     
     q1 = mgs(q1);