[y,Fs] = audioread('mamma_mia.mp3','native');

[phi1,tau]=lorenz1([10 28 8/3 3 -0.1 1 500]);
start=length(y)/2;
max1=length(phi1);
sound(y(start:start+max1),Fs);

pause
scale=0.05;
y1=scale*y(start:start+max1-1)+phi1';
sound(y1, Fs)

pause

x=y1;  %occhio deve essere y1 e non phi1 
f=inline('10*(v-u)','u', 'v','w','x');
g=inline('x*(28-20*w)-v','u','v','w','x');
h=inline('5*x*v-8/3*w','u','v','w','x');

    u(1)   = 1;
    v(1)   = 1; 
    w(1)   = 20;
    
    N      = max1;
    
    
        
                   
                        % Runge-Kutta method: errore globale O(h^4)

                   for i=1:N-1
                           
                          DeltaT=tau(i+1)-tau(i);

                          k1x = DeltaT*f(u(i),v(i),w(i),x(i));
                          k1y = DeltaT*g(u(i),v(i),w(i),x(i));
                          k1z = DeltaT*h(u(i),v(i),w(i),x(i));
                          
                          k2x = DeltaT*f(u(i)+k1x/2, v(i)+k1y/2, w(i)+k1z/2, (x(i)+x(i+1))/2);
                          k2y = DeltaT*g(u(i)+k1x/2, v(i)+k1y/2, w(i)+k1z/2, (x(i)+x(i+1))/2);
                          k2z = DeltaT*h(u(i)+k1x/2, v(i)+k1y/2, w(i)+k1z/2, (x(i)+x(i+1))/2);
                          
                          k3x = DeltaT*f(u(i)+k2x/2, v(i)+k2y/2, w(i)+k2z/2, (x(i)+x(i+1))/2);
                          k3y = DeltaT*g(u(i)+k2x/2, v(i)+k2y/2, w(i)+k2z/2, (x(i)+x(i+1))/2);                     
                          k3z = DeltaT*h(u(i)+k2x/2, v(i)+k2y/2, w(i)+k2z/2, (x(i)+x(i+1))/2); 
                           
                          k4x = DeltaT*f(u(i)+k3x,v(i)+k3y, w(i)+k3z, x(i+1));
                          k4y = DeltaT*g(u(i)+k3x,v(i)+k3y, w(i)+k3z, x(i+1));
                          k4z = DeltaT*h(u(i)+k3x,v(i)+k3y, w(i)+k3z, x(i+1));
                  
                          
                          u(i+1)=u(i)+(k1x+2*k2x+2*k3x+k4x)/6;
                          v(i+1)=v(i)+(k1y+2*k2y+2*k3y+k4y)/6;
                          w(i+1)=w(i)+(k1z+2*k2z+2*k3z+k4z)/6;
                          
                          
             
                          if mod(i,10000)==0
                              i
                          end
                   end    
       

FINALE=i
pause
y2=(y1-u)/scale;
sound(y2, Fs)