///withoutPrompt
function Xd=simul(t,X)
  x=X(:,1);dx=X(:,2:$)
  [xd,dxd]=full_model(x,dx)
  Xd=[xd,dxd];
endfunction

function [xd,dxd]=full_model(x,dx)
  [lhs,rhs]=argn(0);
  z=x(1);th=x(2);zd=x(3);thd=x(4);
  Mati=inv([M+m,m*l*cos(th);cos(th),l]);
  ydd=Mati*[m*l*sin(th)*thd^2+U;g*sin(th)];
  zdd=ydd(1);thdd=ydd(2);
  xd=[zd;thd;zdd;thdd];
  F=[0,m*l*thdd*sin(th)+m*l*thd^2*cos(th);0,zdd*sin(th)+g*cos(th)];
  G=[0,2*m*l*thd*sin(th);0,0];
  if rhs == 2 then dxd=[zeros(2,2),eye(2,2);Mati*F,Mati*G]*dx;end
endfunction ///\withPrompt{}