clear all
f=inline('10*(y-x)','x', 'y','z');
g=inline('x*(28-z)-y','x','y','z');
h=inline('x*y-8/3*z','x','y','z');



store_fin=[0 0 0];

            figure(1)
            a=load('storeLorenz');
            plot3(a(:,1),a(:,2),a(:,3),'b')
            grid on
            hold on
    
            pause
            
for k   = 1:50

    k
    x(1)   = 1+0.001*randn(1);
    y(1)   = 1+0.001*randn(1); 
    z(1)   = 20+0.001*randn(1);
    t      = 0;
    Tfin   = 20;
    DeltaT = 0.05;
    N      = round(Tfin/DeltaT);
    
    
        
                   
                        % Runge-Kutta method: errore globale O(h^4)

                   for i=1:N
                          

                          k1x = DeltaT*f(x(i),y(i),z(i));
                          k1y = DeltaT*g(x(i),y(i),z(i));
                          k1z = DeltaT*h(x(i),y(i),z(i));
                          
                          k2x = DeltaT*f(x(i)+k1x/2, y(i)+k1y/2, z(i)+k1z/2);
                          k2y = DeltaT*g(x(i)+k1x/2, y(i)+k1y/2, z(i)+k1z/2);
                          k2z = DeltaT*h(x(i)+k1x/2, y(i)+k1y/2, z(i)+k1z/2);
                          
                          k3x = DeltaT*f(x(i)+k2x/2, y(i)+k2y/2, z(i)+k2z/2);
                          k3y = DeltaT*g(x(i)+k2x/2, y(i)+k2y/2, z(i)+k2z/2);                     
                          k3z = DeltaT*h(x(i)+k2x/2, y(i)+k2y/2, z(i)+k2z/2); 
                           
                          k4x = DeltaT*f(x(i)+k3x,y(i)+k3y,z(i)+k3z);
                          k4y = DeltaT*g(x(i)+k3x,y(i)+k3y, z(i)+k3z);
                          k4z = DeltaT*h(x(i)+k3x, y(i)+k3y, z(i)+k3z);
                  
                          
                          x(i+1)=x(i)+(k1x+2*k2x+2*k3x+k4x)/6;
                          y(i+1)=y(i)+(k1y+2*k2y+2*k3y+k4y)/6;
                          z(i+1)=z(i)+(k1z+2*k2z+2*k3z+k4z)/6;
                          
                          t=[t i*DeltaT];
    
                   
                          
           
             
                   end    
       
                   
   
            
          %plot3(x,y,z)
            plot3(x(end),y(end),z(end),'r.','MarkerSize',30)

    
           




end