function y=aplusb_cauchy_euler(a,aplusb,c,q1);
% - March 13 2003.(RMM)
% - The aplusb_cauchy_euler(a,aplusb,c,q1) 

%     solves  a x^2 y''+ (a+b) x y' +c y =q1(x)
%      or        a D2 Y +b DY +c Y =q1(exp(t))
% - a , aplusb, and c are constants.  q1 is a function of x
% - The classification is used to find y1 and y2 then 
% - the method of variation is used to evaluate yp.
% - example:
% - usage a=1;aplusb=5;c=-5;q1=x;y=aplusb_cauchy_euler(a,aplusb,c,q1)
% -
   syms x m m1 m2 yh y1 y2 delta c1 c2 alpha beta u1 u2 yp w
  disp(' Second order Linear homogeneous ODE with constant coefficients')
  disp(' a D2Y + b  DY + C Y= q1(exp(t))')
      b=aplusb-a;
sprintf(' where a = %d',a)
    sprintf('       b = %d',b)      
    sprintf(' and,  c = %d',c)

  
      disp('The characteristic equation is  a m^2+bm+c=0')
      %you could use solve here ...    solve('a*m^2+b*m+c', 'm');
         delta = b^2-4*a*c;
      sprintf(' delta= b^2-4ac =%d ',delta)
   
      if delta==0  
         disp(' class B (critical double roots m1=m2=-b/(2a)')
         m1=-b/(2*a);
         m2=m1;
         y1=x^(m1)
         y2=y1*log(x)
      end
      
      if delta >0
            disp(' class A supercritical, two different real roots,')
            m1= (-b+sqrt(delta))/(2*a);
            m2= (-b-sqrt(delta))/(2*a);
            sprintf(' m1=%d',m1)
            sprintf('m2=%d',m2)
            y1=x^m1
            y2=x^m2
         end
         if delta <0
         disp(' class C subcritical , complex conjucate roots')
            alpha=-b/(2*a);
            beta=sqrt(-delta)/(2*a);
            y1=x^alpha*cos(beta*log(x))
            y2=x^alpha*sin(beta*log(x))
         end
                 yh=c1*y1+c2*y2;
         disp('The solutions 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*q1/(a*x^2*w),x)
         u2=int(y1*q1/(a*x^2*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
      return