function out1 = vdp1(t,y)
P = 1;
S0 = pi*0.003^2;
g = 0.01;
S = pi*0.0285^2;
h = 0.13;
m = 0.115;
l = 0.58;
b = 0.1;
lcm = 0.45;
K = 19;
M = 0.350;
out1 = [y(3); y(4); -1/(1/3*m*l^2+(1/12*b^2+l^2)*(P*S0^2*g/(2*S)*t^2-P*S0*sqrt(2*g*h)*t+P*S*h))*(-1*(P*S0^2*g/(2*S)*t^2-P*S0*sqrt(2*g*h)*t+P*S*h)*g*y(1)*lcm+K*h^2*(y(1)-y(2)));
   1/(1/3*m*l^2+1/12*M*b^2+M*l^2)*(-M*g*y(2)*lcm + K*h^2*(y(1)-y(2)))];
