function  y= variation(p,q,r,y1,y2)
% - Revision of March 1st 2002(RMM) 
% - The method of variation of parameters used to solve
%  -  y'' +p y'+q y=r.
% - p , q and r are functions of x.
% - y1 and y2 will be checked if they are true solutions of the homogeneous eqn.
% -   If so then 
% -      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=variation(p,q,r,y1,y2)
syms x c1 c2
    disp(' variation of parameters'); disp(q) 
  disp(' we need to check if y1 and y2 are solutions. ')

            disp(' You are testing  D2y+pDy+q y=0')
            disp(' where p = '); disp(p)
            disp(' and,  q = '); disp(q) 
            %  -  let us also veriy that y1 and y2 are solutions.
             z1=diff(diff(y1,'x'),'x')+p*diff(y1,'x')+q*y1;
             disp(' if z1= 0 then y1 is a solution . z1=');disp(z1)
             z2=diff(diff(y2,'x'),'x')+p*diff(y2,'x')+q*y2;
             disp(' if z2= 0 then y2 is a solution . z2= ');disp(z2)


      if z1==0 
         % - y1 is a solution
      if z2==0
         % - y1 and y2 are solutions
         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=simplify(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(' y2 is not a solution')
end
else
   disp('y1 is not a solution')
   disp('  try again')
end
           return