function yp = lorenz_solver(t,y,flag,r)
        sig = 10;
        b = 8/3;
        yp1= -sig.*(y(1) - y(2));
        yp2= -y(2) +r.*y(1) -y(1).*y(3);
        yp3= -b.*y(3) +y(1).*y(2);
        yp = [yp1;yp2;yp3];
   end