function hamilton
clc;clear all;

[T,Y] = ode15s(@Hamiltonian,[0 5],[0 0 0 0]);

subplot(3,1,1);
plot(T,Y(:,1)*180/pi)
grid;
title('Results of Teta');
ylabel('Teta(1)');

subplot(3,1,2);
plot(T,Y(:,2)*180/pi)
grid;
ylabel('Teta(2)');


subplot(3,1,3);
plot(T,Y(:,1)-Y(:,2))

function dy = Hamiltonian(t,y)
dy = zeros(4,1);
m=1;
g=9.81;
l=1;
n=100;
p=3;
iy=1.46;

dy(1) = (1/(m*l^2*(16-9*cos(y(2)-y(1))^2))) * (12*y(3) - 18*cos(y(2)-y(1))*y(4));
dy(2) = (1/(m*l^2*(16-9*cos(y(2)-y(1))^2))) * (48*y(4) - 18*cos(y(2)-y(1))*y(3));
dy(3) = (-18/(m*l^2*(16-9*cos(y(2)-y(1))^2)^2))*sin(y(2)-y(1))*(6*cos(y(2)-y(1))*y(3)^2 - (16+9*cos(y(2)-y(1))^2)*y(3)*y(4) + 24*cos(y(2)-y(1))*y(4)^2) - 1.5*m*g*l*sin(y(1)) + p*(n*exp(-n*t))*l*cos(y(1));
dy(4) = (+18/(m*l^2*(16-9*cos(y(2)-y(1))^2)^2))*sin(y(2)-y(1))*(6*cos(y(2)-y(1))*y(3)^2 - (16+9*cos(y(2)-y(1))^2)*y(3)*y(4) + 24*cos(y(2)-y(1))*y(4)^2) - 0.5*m*g*l*sin(y(2)) + p*(n*exp(-n*t))*((iy-l)*cos(y(2)));