function y=main()
  
  a=0.001;
  R=1.1;
  c=3;
  m=0.03;
  
  Peq=(log(R)/a)^(1/(1-m))
  Heq=Peq/(c*(1-exp(-a*Peq^(1-m))))
  
  epsi=0.01;
  
  H(1)=(1+epsi)*Heq;
  P(1)=(1+epsi)*Peq;

  
  N=60;

  for n=1:N-1
    H(n+1)=R*H(n)*exp(-a*P(n)^(1-m));
    P(n+1)=c*H(n)*(1-exp(-a*P(n)^(1-m)));
  end
  
  figure(1)
  clf
  hold on
  plot(Heq,Peq,'*k')
  plot(H,P,'*-')
  plot(H(1),P(1),'*-r')
   
  end
  
    