qmax=input(' qmax? ')
dq=qmax/100
q=-qmax:dq:qmax;
pq=size(q);nq=pq(2)
pmax=input(' pmax?')
dp=2*pmax/100
p=-pmax:dp:pmax;
pp=size(p);np=pp(2)

g=1;
V=q.^2*(g/2);
div=(q.^2+1)*2;
coe=1./div;

for i=1:np
for j=1:nq
H(i,j)=p(i)^2*coe(j)+V(j);
end
end
figure(1)
hold off
Z=-2:0.2:+2;
% piano q,p
contour(q,p,H,Z)
xlabel('q')
ylabel('p')
title('Funzione di Hamilton')
% stampa figura
%print -dpsc colham_lev53.ps
%print -dps bwham_lev53.ps

% piano q, qp
qpmax=input('max di qdot? '); 
dqp=2*qpmax/100   
qp=-qpmax:dqp:qpmax;
pqp=size(qp);nqp=pqp(2)

divl=div*(1/4);

for i=1:nqp
for j=1:nq
  E(i,j)=qp(i)^2*divl(j)+V(j);
end
end
figure(2)
hold off
contour(q,qp,E,Z)
xlabel('q')
ylabel('qdot')
title('Funzione energia')
% stampa figura
%print -dpsc ham_lev53b.ps







