function y=cauchy_euler(a,b,c,q1);
% - Revision on March 13 2003.(RMM)
% - The cauchy-euler(a,b,c,q1) solves  a x^2 y''+ b x y' +c y =q1(x)
% - a  b c are constants and 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.
% - 
% - usage a=1;b=3;c=1;q1=x^2;y=cauchy_euler(a,b,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 x^2 D2y + b x Dy + C y= q1')

    sprintf(' where a = %d',a)
    sprintf('       b-a = %d',b-a)      
    sprintf(' and,  c = %d',c)

  
      disp('The characteristic equation is  a m^2+(b-a)m+c=0')
      %you could use solve here ...    solve('a*m^2+(b-a)*m+c', 'm');
         delta = (b-a)^2-4*a*c;
      sprintf(' delta= (b-a)^2-4ac =%d ',delta)
   
      if delta==0  
         disp(' class B (critical double roots m1=m2=-b/(2a)')
         m1=-(b-a)/(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-a)+sqrt(delta))/(2*a);
            m2= (-(b-a)-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-a)/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