function  y= varROO(p,q,r,y1)
% - Revision of March 1st 2002(RMM) 
% - The method of variation of parameters and reduction of order jointly used
%  - to solve y'' +p y'+q y=r.
% - p , q and r are functions of x.
% - y1 is a known solution. 
% -  ROO (reduction of order will find y2 
% -  Then Variation of parameters is used to find y-particular yp. 
% -      yp=y1*u+y2*v
% - where           u=int(-y2*r/w,x)
% -                 v=int(y1*r/w,x)
% -   y=c1 y1 +c2 y2 +yp
% - usage y=(p,q,r,y1)
syms x c1 c2
    disp(' variation plus reducction of order'); disp(q) 
  disp('you know y1 is a solution for ODE . Find y2=y1*u as a second solution. ')

            disp(' You are solving  D2y+pDy+q y=0')
            disp(' where p = '); disp(p)
            disp(' and,  q = '); disp(q) 
            %  -  let us also veriy that y1 is a solution.
            
             z1=diff(diff(y1,'x'),'x')+p*diff(y1,'x')+q*y1;
             disp(' if z1= 0 then y1 is a solution . z1=');disp(z1)
      if z1==0 
       % - y1 is a solution 
         z2=int(p,x)
         z2=exp(-z2)
         z2=simplify(z2/y1/y1)
         u=simplify(int(z2,x))   
         disp(' The second linearly inependent solution y2 is= ')    
         y2=simplify(u*y1)
             yh=c1*y1+c2*y2;
         disp(' the solutiions y1 and y2 are linearly independent.')
         disp(' The wronskian is given by w=y1*diff(y2,x)-y2*diff(y1,x)')
         w=y1*diff(y2,x)-y2*diff(y1,x);
         disp(' the particular solution is given by the method of variation')
         u1=int(-y2*r/w,x)
         u2=int(y1*r/w,x)
         yp = simplify(y1*u1+y2*u2)
         fprintf(' The particular solution')
         
         fprintf(' The Homogeneous solution yh=c1 y1+c2 y2 is ')
         disp(yh)
         
         disp('The  general solution is given by y=yh+yp =')
         y=yh+yp

             
else
         disp(' y1 is not a solution try again')
     end
           return