clear all



f=inline('-x+a*y+x^2*y','x', 'y', 'a','b');
g=inline('b-a*y-x^2*y','x','y','a','b')
a=0.06;
b=0.5;

%err=0;
%h=0;



    
for k   = 1:10 

    k
    x(1)   = 0.5+0.1*k;
    y(1)   = 0.; 
    t      = 0;
    Tfin   = 70;
    DeltaT = 0.1;%10^(-k);
    N      = round(Tfin/DeltaT);

    
        
                   
                        % Runge-Kutta method: errore globale O(h^4)

                   for i=1:N
                          
                          k1x = DeltaT*f(x(i),y(i),a,b);
                          k1y = DeltaT*g(x(i),y(i),a,b);
                          
                          k2x = DeltaT*f(x(i)+k1x/2, y(i)+k1y/2,a,b);
                          k2y = DeltaT*g(x(i)+k1x/2, y(i)+k1y/2,a,b);
                          
                          k3x = DeltaT*f(x(i)+k2x/2, y(i)+k2y/2,a,b);
                          k3y = DeltaT*g(x(i)+k2x/2, y(i)+k2y/2,a,b);                     
                     
                          k4x = DeltaT*f(x(i)+k3x,y(i)+k3y,a,b);
                          k4y = DeltaT*g(x(i)+k3x,y(i)+k3y,a,b);
                          
                  
                          
                          x(i+1)=x(i)+(k1x+2*k2x+2*k3x+k4x)/6;
                          y(i+1)=y(i)+(k1y+2*k2y+2*k3y+k4y)/6;
                          
                          t=[t i*DeltaT];
    
                   
                   
           
             
                   end    
       
                   
            figure(1)
            plot(x,y,'r')
            hold on 
            plot(b,b/(a+b^2),'k*') 
            pause
  %          figure(2)
  %          plot(t,x)
            
%h = [h DeltaT];




end