Heun's method for a ball and beam system

Question

I have a problem with a homework request by the professor at the university but I really can't solve the equation. I have just to resolve the equation using Heun method and choosing appropriately the discretization step (the sampling time interval), adequately motivating the choices made. It's a ball and beam system that after a bit of steps give me these equations:

$\dot x=\begin{bmatrix} \dot x_1\\ \dot x_2\\ \dot x_3 \\ \dot x_4\end{bmatrix}=\begin{bmatrix} x_2\\ x_1x_4^2-g\sin x_3\\ x_4 \\ \frac{-2mx_1x_2x_4-mgx_1\cos x_3+\tau}{mx_1^2+J}\end{bmatrix}= f(x,\tau)$

initial state condition given by the vector:

$ x_0=\begin{bmatrix} x_1 \\ 0 \\ 0 \\ 0\end{bmatrix}$

The only thing I understood is that I have to reduce that to a first grade equation imposing $ \frac{\partial x}{\partial t} = y$ but I don't know how to resolve then. Don't know how to proceed solving that problem with Heun method. Because I never used it and even if I see a lot of video online about it I can't resolve that problem up here. Only started using Matlab but don't know how to use it even if the prof said to us to use it if we want to resolve that. I never used finite difference methods unfortunately...

Answer 1

You write the right side of your equation in Matlab format

  function dotx = ballbeam(t,x)
      dx2 = x(1)*x(4)^2 - g*sin(x(3));
      dx4 = (-2*m*x(1)*x(2)*x(4)-m*g*x(1)*cos(x(3))+tau)/(m*x(1)^2+J);
      dotx = [ x(2), dx2, x(4), dx4 ];
  end

and then feed this into a generic Heun loop

  function x = odeheun(f,t,x0)
      x = [x0];
      for i=1:(length(t)-1)
          h = t(i+1)-t(i);
          k1 = h*f(t(i),x(i,:));
          k2 = h*f(t(i+1),x(i,:)+k1);
          x(i+1,:) = x(i,:)+0.5*(k1+k2);
      end%for
  end%function

and then call it like (tested in octave)

  x0 = [x1_ini, 0, 0, 0 ];
  t = t0:h:tf;
  x = odeheun(@(t,x)ballbeam(t,x),t,x0); % possibly only @ballbeam in matlab
  plot(t,x(:,1));