% Volterra (Ex.5 Kuznetsov 2nd ed., p. 32)
  close all
  clear all

% Parameters
  gamma = 1 ; dt = 5e-2 ; nstep = 500 ;

% Grid of initial conditions
  for x0 = [0:1.5:3]
   for y0 =  [0:1.5:3]

    t = 0 ; x = [x0 ; y0] ; xvec = [x0] ; yvec = [y0] ; tvec = [t] ; 
  
    for tstep = 1:nstep
%    Forward step via RK4
      p1 = voltfun(x,gamma) ;
      p2 = voltfun(x+0.5*dt*p1,gamma) ;
      p3 = voltfun(x+0.5*dt*p2,gamma) ;
      p4 = voltfun(x+dt*p3,gamma) ;
  
      x = x + (dt/6)*(p1+2*p2+2*p3+p4) ;
      t = dt*tstep ;
      xvec = [xvec x(1)] ; yvec = [yvec x(2)] ; tvec = [tvec t];
    end

  
     plot(xvec,yvec,'-k','linewidth',1.25) ; drawnow ; hold on
     plot(xvec(1),yvec(1),'ok','markerfacecolor',[1 1 1]) 

end
end


 axis equal % preserve aspect ratio
 axis([0 4 0 4]) ; 
   
%  Isocurves of H(x,y) = ln(y*x^gamma) - (x+y)
   xx = [.25:.01:4.25] ; yy = [.25:.01:4.25] ; 
   [XX,YY] = meshgrid(xx,yy) ;
    ZZ = log((XX.^gamma).*YY)-(XX + YY)  ;
    zmax = max(max(ZZ)) ; zmin = min(min(ZZ)) ; 
    c1 = zmin + (zmax - zmin)*.75 ; c2 = zmax - (zmax -zmin)*.75 ;
    c = linspace(zmin,zmax,7);
    contour(XX,YY,ZZ,c,'color',[.5 .5 .5])

  

