function en221hw4


       B0 = 1.;
       muonrho = 0.4;
       a_initial = 1;
       b_initial = 5.;
       dadt_initial = 0;
       initial_w = [a_initial,dadt_initial,b_initial];
       options = odeset('Events',@event);      
     
       muonrho = 0.1;
       [times,sols] = ode45(@eom,[0,80.2],initial_w,options);
       plot(times,sols(:,1));

       hold on;
     
       muonrho = 0.2;
       [times,sols] = ode45(@eom,[0,80.2],initial_w,options);
       plot(times,sols(:,1));       

       muonrho = 0.4;
       [times,sols] = ode45(@eom,[0,80.2],initial_w,options);
       plot(times,sols(:,1));          
       
       muonrho = 1.;
       [times,sols] = ode45(@eom,[0,80.2],initial_w,options);
       plot(times,sols(:,1));   
       
       muonrho = 2.;
       [times,sols] = ode45(@eom,[0,80.2],initial_w,options);
       plot(times,sols(:,1));  
       
    function dwdt = eom(t,w)
       a = w(1); dadt = w(2); b = w(3); 
       
       d2adt2 = -B0*b/2/a - 4*dadt*muonrho*(a^2+b^2+a*b)/a^2/b^2 ...
       -dadt^2*(b-a)*(a^2+2*a*b+3*b^2)/2/a/b^3;
       
       dbdt = dadt*a^2/b^2;
       
       dwdt = [dadt;d2adt2;dbdt];
   
    end

    function [ev,stop,dir] = event(t,w)
        ev = w(1)-0.01;
        dir = 0;
        stop = 1;
    end


end

