Runge-Kutta Method solving four ODES

Question

I am currently working on some system of ODEs enter image description here.

These four odes are equation of motion for a physical system, photo linked below:

Scroll down to the bottom of this page, the physical system in motion and interactive java applet: https://www.myphysicslab.com/pendulum/cart-pendulum-en.html.

x represents displacement of the main block
q represents the angle theta of the pendulum from the vertical position
v represents velocity of the main block
w represents the angular velocity omega of the pendulum.

I coded myself a matlab file to draw a graph of x vs t(time), however, it did not work and produced extremely large outputs. I wonder if anyone could help me with the code (since I am a newbie) and give me some advices.

Btw, I chose to use this numerical method because the equations give implicit solutions when solved normally (nonlinear).

% x'= v
% q'= w
% v' = f(x,q,v,w) = (p*R*(w^2)*sin(q)+p*g*sin(q)*cos(q)-k*x-c*v+(b/R)*w*cos(q))/((m+p*(sin(q))^2))
% w' = g(x,q,v,w) = (-p*R*(w^2)*sin(q)*cos(q)-(p+m)*g*sin(q)+k*x*cos(q)+c*v*cos(q)-(1+M/q)*(b/R)*w)/(R*(M+p*(sin(q))^2))

p = 0.05;     % mass of pendulum mass
m = 0.198;    % mass of large mass
R = 0.1;      % length of pendulum rod
g = 9.80665;  % gravitational acceleration
k = 15;       % large mass spring constant
c = 0.124;    % large mass damping coefficient
b = 0.00543   % damping torque coefficient for pendulum mass

% Equations
fX = @(t,V) V;
fQ = @(t,W) W;
fV = @(t,X,Q,V,W) (p*R*(W^2)*sin(Q)+p*g*sin(Q)*cos(Q)-k*X-c*V+(b/R)*W*cos(Q))/((m+p*(sin(Q))^2));
fW = @(t,X,Q,V,W) (-p*R*(W^2)*sin(Q)*cos(Q)-(p+m)*g*sin(Q)+k*X*cos(Q)+c*V*cos(Q)-(1+m/p)*(b/R)*W)/(R*(m+p*(sin(Q))^2));

% Initial conditions
t(1) = 0;
X(1) = 0.1;
Q(1) = 0;
V(1) = 0;
W(1) = 0;

% Step size and number of steps
h = 0.0001;
time_total = 20;
n = ceil(time_total/h);

% RK4 loop
for i = 1:n    
    t(i+1) = t(i) + h; % time increment
    
    %RK1
    k1X = fX(t(i),V(i));
    k1Q = fQ(t(i),W(i));
    k1V = fV(t(i),X(i),Q(i),V(i),W(i));
    k1W = fW(t(i),X(i),Q(i),V(i),W(i));
    
    %RK1
    k2X = fX(t(i)+h/2,V(i)+h/2+k1V);
    k2Q = fQ(t(i)+h/2,W(i)+h/2+k1W);
    k2V = fV(t(i)+h/2,X(i)+h/2+k1X,Q(i)+h/2+k1Q,V(i)+h/2+k1V,W(i)+h/2+k1W);
    k2W = fW(t(i)+h/2,X(i)+h/2+k1X,Q(i)+h/2+k1Q,V(i)+h/2+k1V,W(i)+h/2+k1W);
    
    %RK1
    k3X = fX(t(i)+h/2,V(i)+h/2+k2V);
    k3Q = fQ(t(i)+h/2,W(i)+h/2+k2W);
    k3V = fV(t(i)+h/2,X(i)+h/2+k2X,Q(i)+h/2+k2Q,V(i)+h/2+k2V,W(i)+h/2+k2W);
    k3W = fW(t(i)+h/2,X(i)+h/2+k2X,Q(i)+h/2+k2Q,V(i)+h/2+k2V,W(i)+h/2+k2W);
    
    %RK1
    k4X = fX(t(i)+h,V(i)+h+k3V);
    k4Q = fQ(t(i)+h,W(i)+h+k3W);
    k4V = fV(t(i)+h,X(i)+h+k3X,Q(i)+h+k3Q,V(i)+h+k3V,W(i)+h+k3W);
    k4W = fW(t(i)+h,X(i)+h+k3X,Q(i)+h+k3Q,V(i)+h+k3V,W(i)+h+k3W);
    
    %Summing
    X(i+1) = X(i)+(h/6)*(k1X+2*k2X+2*k3X+k4X);
    Q(i+1) = Q(i)+(h/6)*(k1Q+2*k2Q+2*k3Q+k4Q);
    V(i+1) = V(i)+(h/6)*(k1V+2*k2V+2*k3V+k4V);
    W(i+1) = W(i)+(h/6)*(k1W+2*k2W+2*k3W+k4W);
    
end

plot(t,X);

Answer 1

One obvious problem is that

k2X = fX(t(i)+h/2,V(i)+h/2+k1V);

should read

k2X = fX(t(i)+h/2,V(i)+h/2*k1V);

multiplication, not addition of $h/2$.

The general critique is that you are using the Matrix Laboratory, so in keeping with its philosophy, use vector valued functions and the built-in vector operations. Also use the provided ODE integrators and make your interfaces similar to them so that switching back and forth is easy.