function dv = dof3(t,v)
dv = zeros(4,1); 
m   =3.085; % mass, tonnes
m_f =0.060; % mass of front wheel, tonnes
m_r =0.060; % mass of rear wheel, tonnes
I_xxb=1.190; % sprung mass roll inertia, tm^2
I_zzb=7.050; % sprung mass yaw inertia, tm^2
l=3.650; % a+b wheelbase, m
b=1.520; % b CoG distance from rear axle, m
a=l-b; % a 
Cf=140; % Front cornering stiffness, kN/rad
Cr=190; % Rear cornering stiffness, kN/rad
U=16.66; % Forward speed,m/s
I_zzr=0.0012; % yaw inertia of rear wheel, tm^2
I_zzf=0.0012; % yaw inertia of front wheel, tm^2
mb = m-2*0.160 ; % sprung body mass, 
h_cf= 0.760  ; % Height of roll centre, m
h_cr= 0.760  ; % Height of roll centre, m
h_a=  0.260   ; % height of reference roll axis, m
K_phi=10000  ; % Torsional stiffness coefficient, kNm/rad
C_phi=0.250  ; % Torsional damping coefficient, 
h=0.5 ;
D_phi=0.5*C_phi*0.1*0.1 ; 
I_xx = I_xxb + mb*(h)^2;
I_zz = I_zzb + I_zzf + I_zzr + m_f*a^2 + m_r*b^2;
df = h_cf-h_a;
dr = h_cr-h_a;
g=9.8;
ass=0.3*9.81; % Steady state acceleration, m/s^2 
delf=(l^2*Cf*Cr+m*U^2*(b*Cr-a*Cf))*ass/(l*Cf*Cr*U^2); % Steer angle front, rad

Mass = [m a*m_f-b*m_r mb*h 0; a*m_f-b*m_r I_zz 0 0;mb*h  0 I_xx 0; 0 0 0 1];
M1   =-[(Cf+Cr)/U (a*Cf-b*Cr)/U+m*U 0 0;(a*Cf-b*Cr)/U ((a^2)*Cf+(b^2)*Cr)/U+(a*m_f-b*m_r)*U 0 0;(df*Cf+dr*Cr)/U (a*df*Cf-b*dr*Cr)/U+mb*h*U D_phi K_phi-mb*g*h;0 0 -1 0];
M2   = [Cf*delf;a*Cf*delf;df*Cf*delf;0];

dv = inv(Mass)*M1*v+inv(Mass)*M2 ;  % equation

end