Pendulum on the spring in 3D. 
Dynamic simulation model.
The position of the hanging is externally controled

The externally controled hanging position has coordinates {xp, yp, zp}
The spring has length at rest l0 and spring constant K
The pendulum mass is m

The parametrization uses Cardan like parametrization with first rotation around x axis by angle α and then rotation around new z axis by the angle β.
This parametrization is used to avoid singularity in vertical position of the pendulum.
In fact this parametrization has another singularity when pendulum is in y axis direction.

Independent variables

vars = {l[t],α[t],β[t]};

Bounding equations

x=l[t] Cos[β[t]]Cos[α[t]] +xp[t];
y = l [t]Sin[β[t]] +yp[t];
z=-l [t]Cos[β[t]]Sin[α[t]] +zp[t];

First and second derivatives of bounding equations

D[x,t]//TraditionalForm
D[y,t]//TraditionalForm
D[z,t]//TraditionalForm
D[x,{t,2}]//Expand//TraditionalForm
D[y,{t,2}]//Expand//TraditionalForm
D[z,{t,2}]//Expand//TraditionalForm

Force components

Fx[t]:=Force[t] Cos[β[t]]Cos[α[t]]
Fy[t]:=Force[t] Sin[β[t]]
Fz[t]:=Force[t] Cos[β[t]]Sin[α[t]]
Force[t]:=K (l[t]-l0)

Newton equations

eq1=m D[x,{t,2}]+Fx[t];
eq2 = m D[y,{t,2}] +Fy[t];
eq3 = m D[z,{t,2}] +G-Fz[t];
eq = {eq1,eq2,eq3};

Left sides of  Newtonov equations in readable format

eq1 //Expand//TraditionalForm
eq2 //Expand//TraditionalForm
eq3 //Expand//TraditionalForm

Newton equations could be rewritten into the form A w = b and numerically solved: w = (A^(-1)) b

Vector of second derivatives of independent variables w

w = D[vars,{t,2}]

Matrix A (using MapAt)

A = (Coefficient[eq,#]&/@w )\[Transpose];
A//TraditionalForm

Vector of right sides b

b=A . w - eq//Simplify;
b //Expand//TraditionalForm

We could make inversion of such a simple system symbolically but standard procedure would be to make inversion numerically in Matlab.

AI =Inverse[A]//FullSimplify;
AI //TraditionalForm

AI.b//FullSimplify //TraditionalForm

Preparations for solutions in Matlab

Finally we convert a system of 3 partial diferencial second order equations to the system of 6 partial diferencial first order equation using subsitution e.g.
z1 = l
z2 = α
z3 = β
z4 = l’
z5 = α’
z6 = β'
which allows to use function ode45 in  Matlab to solve equations while knowing initial values w0 and external input {xp[t], yp[t], zp[t]}
The equations will be entered into ode45 as

Final equations for ode45

eqfin = Flatten[{l',α',β',AI.b}]//FullSimplify //TraditionalForm

Conversion to matlab notation

http : // library.wolfram.com/infocenter/MathSource/577/
modify the path at the next line

Get["ToMatlab.m",Path->"skola/Math/."]

Possible name conversions, the next expression says, how to name variables in Matlab. 
  One can easily modify the names according one' s choice.
  Also is possible to leave t as a parameter of the functions e.g. by rules {α -> alpha, ...}, 
which will result in "alpha(t)".

cond={l[t]->l,α[t]->alpha,β[t]->beta,l'[t]->dl,α'[t]-> dalpha,β'[t]->dbeta,
xp''[t]->ddxp,yp''[t]->ddyp,zp''[t]->ddzp};

ToMatlab looks similar but leave unwanted hidden characters '\'.
    PrintMatlab works smoothly when you want copy - paste the result into a Matlab.

PrintMatlab[A/.cond]
PrintMatlab[b/.cond]