%  function f=rod(t,x,y)
function f=rod(t,X) 
% Лично мне привычнее сводить все в вектор состояния 
% Тогда система в матричном виде будет X'=A(X,t,...)*X 
%     X=[x
%           x'
%           y
%           y']
x=X(1:2);
y=X(3:4);
    lambda=2*10^7;
    l=100;
    d=6;
    theta=1;
    F=30;
    E=2*10^6;
    k=l^2/(d*pi^2);
    h=1;%  Я решил, что можно взять за 1
    
    f=[x(2)
        (4*E*F*lambda*theta-8*lambda^2*k*t^2+...
        2*E*F*h*lambda+4*E*F*k*lambda*y(1)+...
        4*E*F*k*lambda*x(1)+8*E*F*k*lambda*t*y(2)+...
        8*E*F*k*lambda*t*x(2))/(E^2*F^2*h-...
        4*E*F*k*lambda*t^2)
       y(2)
        (8*k*lambda^2*t^2+4*E*F*lambda*theta-...
        2*E*F*h*lambda+4*E*F*k*lambda*y(1)+...
        4*E*F*k*lambda*x(1)+8*E*F*k*lambda*t*y(2)+...
        8*E*F*k*lambda*t*x(2))/(E^2*F^2*h-...
        4*E*F*k*lambda*t^2)];