
(* :Name: Calculus`DSolve` *)

(* :Title: Differential Equations of the Hypergeometric Type and
           Selected Nonlinear firsti/second order equations *)

(* :Author: Victor S. Adamchik and Alexei V. Bocharov respectively*)

(* :Context: System` *)

(* :Package Version: 2.1 *)

(* :Copyright: Copyright 1992  Wolfram Research, Inc.

	Permission is hereby granted to modify and/or make copies of
	this file for any purpose other than direct profit, or as part
	of a commercial product, provided this copyright notice is left
	intact.  Sale, other than for the cost of media, is prohibited.

	Permission is hereby granted to reproduce part or all of
	this file, provided that the source is acknowledged.
*)

(* :History:
	Version 2.1 by Victor S. Adamchik,  1991 and
	            by Alexei V. Bocharov,  1992.
                    Last touch by A.V.Bocharov: March 1992.
*)

(* :Keywords:  *)

(* :Source:
	V. Adamchik, The Dissertation Thesis, 1987,
	A. Bocharov, Unpublished
*)

(* :Mathematica Version: 2.1 *)

(*****************************************************************************)
BeginPackage["Calculus`DSolve`"]

Unprotect[ DSolve ]

DSolve::dsing =
"Unable to fit initial condition `1` at singular point `2`."

DSolve::dvsep =
"Unable to separate variables due to complexity of intermediate expressions."

Begin["`Private`"]

(*****************************************************************************)

DSolve[prob_,uu_,t_] :=
      Module[ {DS,memo,DSFail},


DS[y_'[x_]==F_,y_[x_]] :=

     Module[

	    { wF , var, key , ShiftFactor, SepVarQ , CautionSolve,
              IMultiplier, ScaleFactor, DeepThought, Reducible1,
	      Reducible2, GeneralizedBernoulli, GeneralizedScale,
	      DeeperThought}

	    , (* here goes the second argument of the Module *)

	    ShiftFactor :=
	      Module [ {w} , w = Simplify[ D[wF,x] / D[wF,var] ] ;
		       If [ FreeQ[w,var] && FreeQ[w,x] ,
			   w, Infinity ] ] (*module*) ;

            ScaleFactor :=
	       Module [ {w},
			w = Simplify[Expand[
			( wF - var D[wF,var]) / (wF + x D[wF,x])]];
			If [ FreeQ[w,var] && FreeQ[w,x] ,
			     w,
			     Infinity ]
                      ] ; (*eludom*)

            DilatationFactor :=
	       Module [ {w, Den },
			Den = Simplify[Expand[x D[wF,var] -1 ]] ;
			If [ Den === 0 ,
			     Infinity ,
			      w = Simplify[Expand[
				   ( wF^2 + var D[wF,x])/ Den ]];
                              If [ FreeQ[w,var] && FreeQ[w,x] ,
				   w,
				   Infinity
                                 ]
                            ]
                       ] ; (*eludom*)

            ScaleFactors :=
               Module[ { a , b , c , Den, Fx, Fy, Fxx, Fxy, Fyy },
			  Fx=D[wF,x];
			  Fy=D[wF,var];
			  Fxx=D[Fx,x];
			  Fxy=D[Fx,var];
			  Fyy=D[Fy,var];
			  Den = Simplify[ -(Fx*Fxy*Fy) + Fx^2*Fyy +
			  Fxy^2*wF - Fxx*Fyy*wF] ;
                          If[ Den === 0,
			     Infinity,
			        (* Comments: zero Den has a very complicated  *)
				(*           meaning.                         *)
				(* If we note, however that Den = Fx D1 + F D2*)
				(* where                                      *)
				(*      D1 = Fx Fyy - Fy Fxy ,                *)
				(*      D2 = Fxx Fyy - Fxy^2 ,                *)
                                (* Then we can go into the special case       *)
				(*      D1 === 0 , D2 === 0                   *)
				(* the contents of which is that either       *)
				(* F === G[ y + a x + b]                      *)
				(* and thus a shiftable equation              *)                                (* OR                                         *)                                (*    F is linear in y[x]                     *)
				(*                                            *)                                (* ****************************************** *)

                             a = Simplify[
			     (Fx*Fy^2 - Fxy*Fy*wF + Fx*Fyy*wF + Fx*Fxy*Fy*x -
			     Fx^2*Fyy*x - Fxy^2*wF*x + Fxx*Fyy*wF*x)/Den];
			     If[ FreeQ[a,x] && FreeQ[a,var],
				 b = Simplify[
                                 (-3*Fx*Fxy*Fy + Fxx*Fy^2 + 2*Fx^2*Fyy +
				 Fxy^2*wF - Fxx*Fyy*wF)/Den];
				 If[ FreeQ[b,x] && FreeQ[b,var],
				    c = Simplify[(-(Fx^2*Fy) - Fx*Fxy*wF +
				    Fxx*Fy*wF + 3*Fx*Fxy*Fy*var -
				    Fxx*Fy^2*var - 2*Fx^2*Fyy*var -
				    Fxy^2*wF*var + Fxx*Fyy*wF*var )/Den];
				    If[ FreeQ[c,x] && FreeQ[c,var] ,
					{a,b,c},
					(*else*)
					Infinity
                                      ] (*fi*)
				      ,
                                    (*else*)
				    Infinity
                                   ] (*fi*)
				   ,
                                 (*else*)
				 Infinity
                               ] (*fi*)
                        ] (*fi*)
                 ] (*eludom*)

                 ;

            Reducible1[ wF_,var_,z_ ] :=
               Module[ { Fy , Fyy, Den ,f1 },
		       Fy = D[ wF , var ];
		       Fyy = D [ Fy , var ] ;
		       Den = Simplify[Expand[Fyy^2 - Fy D[Fyy, var ] ] ] ;
		       If [ Den === 0 ,
			    Infinity,
			    (*else*)
			    f1 = Simplify[Expand[(-D[Fy,z] Fyy +
					     Fy D[Fyy,z] ) / Den ] ] ;
                            If [ FreeQ[ f1 , var ] ,
				 If [
				   Simplify[ Expand[
				      D[f1,z] Fy - f1 Fy^2 +
				      wF f1 Fyy - f1^2 Fyy -
				      Fy D[wF,z] + (wF - f1) D[Fy,z]
				      ]] === 0 ,
				      f1,
				      Infinity
                                     ]
                                  ,
				  (*else*)
				  Infinity
                                ] (*fi*)
                            ] (*fi*)
                      ] (*eludom*)  ;

            Reducible2[ wF_,var_,z_ ] :=
	       Module[ { u, g , Den , Fy, Fyy , Fyyy, Fx, Fxy , Fxyy },
		       Fy = D[wF, var ] ;
		       Fx = D[wF, z] ;
		       Fxy = D[Fx , var ] ;
		       Fyy = D[Fy , var ] ;
		       Fxyy = D[Fxy , var ] ;
		       Fyyy = D[Fyy , var ] ;

		       Den = Simplify[ Expand[ 2*wF*Fyy - 2*var*Fy*Fyy +
					  var^2*Fyy^2 + var*wF*Fyyy -
					  var^2*Fy*Fyyy ] ] ;
                       If [ Den === 0 ,
			    Infinity ,   (*This is however a palliative*)
					 (* Den===0 has a meaning !    *)
                            g= Together[ Simplify [ Expand[
				  (Fyy*Fx - var*Fyy*Fxy -
				   wF*Fxyy + var*Fy*Fxyy)/Den ] ] ] ;
                            If [ FreeQ[ g , var] ,
				(*then*)
				g = Simplify[ PowerExpand[Exp[
				       Integrate[ g , z ] ] ] ] ;
                                If [ SepVarQ [ (wF/.var:>(g u )) -
					D[g,z] u , u,z  ] ,
                                     (*then*)
				     g,
				     Infinity
                                   ],
                                (*else*)
				Infinity
                               ]
                           ]
                      ] (*eludom*) ;

            GeneralizedBernoulli[wF_,var_,z_] :=
               Module[ { Fx , Fy , Fxy ,Den, k  } ,
		       Fx= D[wF,z] ;
		       Fy= D[wF,var] ;
		       Fxy= D[Fx,var] ;
		       Den = Fxy wF - Fx Fy ;
                                         (* Recall that Den = 0 means  *)
					 (* equation with separating   *)
					 (* variables : so if separation
					    of variables shadows Bernoulli
					    then the Den must be nonzero
					    here.                      *)
                                         (* SEP VAR MUST SHADOW        *)
					 (* GEN BERNOULLI              *)
                       k=Simplify[(wF D[Fxy,var] - Fx D[Fy,var])/Den];

		       If [ FreeQ[k,z] &&
			    FreeQ[Simplify[Expand[Fy-k wF]],var] ,
			    (*then*)
			    k,
			    Infinity
                          ] (*if*)
                     ] (*eludom*)

                     ;

	    RiccattiSolution[ a_ , b_ , z_ , s_ , var_ ] :=
				         (* s = s[x] is a special solution
					    of the Riccatti equation   *)
                                         (* y'[x] == a y[x]^2 + b y[x] + c
					  *)
	       Module[ { S } ,
		  (* z'[x] == - ( 2 a s + b ) z[x] - a *)
		       S = Integrate[ 2 a s + b , z ];
		       {{var-> s + 1/
		       Simplify[Expand[
			  (C[1] - Integrate[a Exp[S],z]) Exp[-S] ]]
			  }}
                     ] (*eludom*)
		     ;


RiccattiSolve[ wF_ , var_ , z_ ] :=
   Module[ { a , b , c , A , RR , L } ,
	   L = CoefficientList[ Collect[wF,var] ,var];
	   a=Last[L];
	   b=L[[2]];
	   c=First[L];
	   If [ FreeQ[ Simplify[b/a] , z ] &&
		FreeQ[ Simplify[c/a] , z],
		(*then*)
		RiccattiSolution[ a, b, z ,
		 Simplify[
		    PowerExpand[
		       Simplify[( -b + Sqrt[ b^2 - 4 a c ] ) / ( 2 a ) ]]] ,
				  var ] , (* then clause *)
                (*else*)
		RR = Sqrt[ (b z - 1)^2 - 4 a z^2 c ] / ( 2 a z^2 ) ;
		A = Simplify[
		       PowerExpand[
			  Simplify[(1- b z ) / ( 2 a z^2) + RR ]]];
                If [ FreeQ[A,z],
		     RiccattiSolution[ a, b, z, A z , var ]
                     ,
		     (*else*)
		     A = Simplify[
			    PowerExpand[
			       Simplify[(1- b z ) / ( 2 a z^2) - RR ]]];
                     If [ FreeQ[A,z],
			  RiccattiSolution[ a, b, z, A z , var ]
                          ,
			  (*else*)
			  RR = Sqrt[ (b z + 1)^2 - 4 a z^2 c ] / ( 2 a ) ;
			  A = Simplify[
				 PowerExpand[
				    Simplify[(-1- b z ) / ( 2 a ) + RR ]]];
                          If [ FreeQ[A,z],
			       RiccattiSolution[ a, b, z, A / z , var ]
			       ,
			       (*else*)
			       A = Simplify[
				      PowerExpand[
					 Simplify[(-1- b z ) / ( 2 a ) - RR ]]];                               If [ FreeQ[A,z] ,
                                    RiccattiSolution[ a, b, z, A / z , var ]
				    ,
				    (*else*)
				    Infinity
                                  ] (*fi*)
                             ] (*fi*)
                        ] (*fi*)
                   ] (*fi*)
              ] (*fi*)
         ] (*eludom*)
	 ;

	    SepVarQ[wF_,var_,z_] :=
	       Module [ {w} , w = Factor[wF,Trig->True];
			Simplify[D[D[w,z],var] w - D[w,z] D[w,var]] === 0 ]
                       (*eludom*) ;

            CautionSolve [ eq_, term_ ] :=
	       Module [ {ws} ,
		    Off[Solve::ifun];
		    ws = Solve[ eq/.
		    { c_. Log[a_]+d_. Log[b_] :> Log [ a^c b^d ] }, term ] ;
		    On[Solve::ifun];
		    If [ MatchQ[ws,{}],
		    DSFail,
		    ws ] ]
		    (*eludom*) ;

	    IMultiplier [ m_ , wF_ , var_ , z_ ] :=
                                         (* Integrate using Integrating *)
					 (*  Multiplier  m_             *)
               Module [ {  w } ,
		     w= Integrate [ m , var ] ;
		     Simplify[(*Expand*) ( w - Integrate [
			D[w, z ] + wF m , z ] )] == C[1]
                   ] (*eludom*) ;

	    GeneralizedScale[B_] :=
	       Module[ { Fx , Fy , Fxy ,  Den, a, b , wB } ,
		       Fx= D[wF,x] ;
		       Fy= D[wF,var] ;
		       Fxy= D[Fx,var] ;
		       Den = Fxy wF - Fx Fy ;
					 (* Recall that Den = 0 means   *)
					 (* equation with separating    *)
					 (* variables : so if separation*)
					 (* of variables shadows GenScal*)
					 (* the Den must be nonzero here*)

                       b= 1/D[B,var];    (* this is the symmetry coeff. *)
                       a = Simplify[Expand[
			      ( wF^2 D[D[b,var],var] - wF Fy D[b,var] +
			      ( Fy^2 - wF D[Fy,var] ) b ) / Den ]] ;
					 (* this is, potentially, the other *)
                       If[ FreeQ[a,var] ,
			  (*then*)
			  If[Simplify[Expand[
				b Fy + a Fx - D[b,var] wF + D[a,x] wF
				]] === 0 ,
                             (*then*)
			     a,
			     (*else*)
			     Infinity
			     ] (*fi*)
			     ,
                          (*else*)
			  Infinity
                         ] (*fi*)
                      ] (*eludom*)
		      ;

	    DeeperThought :=
					 (* Still more sophisticated    *)
					 (* tricks implemented here     *)
            Module[{a},
			 a=GeneralizedScale[-1/var];
			 If[a=!=Infinity,
			    CautionSolve[
			       IMultiplier[var^2-a wF,wF,var,x],var]/.
			       var:>y[x] ,
                            (*else*)
			    a=GeneralizedScale[-1/var^2];
			    If[a=!=Infinity,
			       CautionSolve[
				  IMultiplier[var^3/2-a wF,wF,var,x],var]/.
				  var:>y[x] ,
                               (*else*)
			       a=GeneralizedScale[-1/var^3];
			       If[a=!=Infinity,
				  CautionSolve[
				     IMultiplier[var^4/3-a wF,wF,var,x],var]/.
				     var:>y[x] ,
                                  (*else*)
				  a=GeneralizedScale[var^2/2];
				  If[a=!=Infinity,
				     CautionSolve[
					IMultiplier[1/var-a wF,wF,var,x],
					var]/.
					var:>y[x] ,
                                     (*else*)
				     a=GeneralizedScale[var^3/3];
				     If[a=!=Infinity,
					CautionSolve[
					   IMultiplier[1/var^2-a wF,wF,var,x],
					   var]/.
					   var:>y[x] ,
                                        (*else*)
				        Infinity
                                       ] (*fi*)
                                    ] (*fi*)
                                 ] (*fi*)
                              ] (*fi*)
                           ] (*fi*)
                      (* ] fi*)
                     (* ] fi*)
		     ] (*eludom*)
                     ;

            DeepThought[wF_ , var_ , z_ ]      :=
					 (* More sophisticated tricks   *)
					 (* implemented here            *)
            Module[ { u , g , uF , w }
                  ,
		  key=Reducible1[wF,var,z];
		  If[key=!=Infinity,
		    (*then*)
                    uF = Simplify[Expand[
			    ( wF/.var:>(Integrate[key,z] + u )) - key
			    ]] ;
                    w=Simplify[ D[uF,z]/uF ] ;
                    g = Exp[Integrate[ w , z ]] ;
		    If [ FreeQ [ g , var ] ,
                        CautionSolve[
			   (Simplify[Integrate[ Simplify[Expand[g/uF]],
					u]/.u:>(var - Integrate[key,z])])
                           == (Integrate[g, z] + C[1]) , var ] (*cauti*)
                        ,
			DSFail
                       ] (*fi*)
                    ,
		    (*else*)
		    key=Reducible2[wF,var,z];
		    If[key=!=Infinity,
		      (*then*)
                      uF = Simplify[Expand[
			      ((wF/.var:>(key u)) - D[key,z] u)/key ]];
                      w=Simplify[ D[uF,z]/uF ] ;
		      g = Exp[Integrate[ w , z ]] ;
		      If [ FreeQ [ g , u ],
			  (*then*)
			  CautionSolve[
			     (Simplify[Integrate[ Simplify[Expand[g/uF]],
				       u]/.u:>(var - Integrate[key,z])])
                             == (Integrate[g, z] + C[1]) , var ] (*cauti*)
                          ,
			  DSFail
                         ] (*fi*)
                       ,
		       (*else*)
		       key=GeneralizedBernoulli[wF,var,z];
		       If[key=!=Infinity,
			 (*then*)
			 g = Simplify[D[wF,var] - key wF];
                         CautionSolve[
			    IMultiplier [
			       1/Simplify[ PowerExpand [ Exp[
				  Integrate[key,var]] *
				  Exp[Integrate[g,z] ]]],wF,var,z],
				  var ] (*cautionsolve*)
                         ,
			 (*else*)
                         If[Simplify[D[D[D[wF,var],var],var]]===0,
			    (*then*)
			    RiccattiSolve[wF,var,z]
			    ,
			    (*else*)
			    Infinity
                           ]
			 ] (*fi*)(*generalized bernoulli clause *)
                     ] (*fi*)(*Reducible2 clause*)
                 ] (*fi*)(*Reducible1 clause*)
               ] (*eludom*) (*DeepThought*)
               ;

(************************** front-end block *****************************)

	    wF=F/.y[x]:>var ;            (* Good y[x] is a dead y[x]    *)

            If[ FreeQ[wF,var] , {{ y[x] -> Integrate[ F , x ] + C[1] }} ,
	    (*else1*)
	      If[ SepVarQ[wF,var,x] ,    (* Separating variables clause *)
                  (*then*)
		  If [ PolynomialQ[wF,var],
		       (*then*)
		       key = Last[CoefficientList[Collect[wF,var],var]];
		       CautionSolve[ Integrate[key/F , y[x] ]
			  == Integrate[key,x ] + C[1],y[x]],
                       (*else3*)
                      If [ PolynomialQ[wF,x],
			   (*then*)
                           key = Last[CoefficientList[Collect[F,x],x]];
                           CautionSolve[ Integrate[ 1/key, y[x] ]
			      == Integrate[ F/key , x ] + C[1],y[x]],
                           (*else4*)
			   key = Exp[Integrate[Simplify[ D[wF,x]/wF ],
				    x ]] ;
                           If[ FreeQ[key,var] ,
                               CautionSolve[ Integrate[Simplify[key/
			       F], y[x] ] == Integrate[key, x] + C[1] ,
			       y[x] ] , (*then CautionSolve*)
			       (*else5*)
			       Message[DSolve::dvsep];
			       DSFail
                             ] (*fi5*)
                         ] (*fi4*)
                     ] (*fi3*)
                  , (* sep var clause ended here *)
		  (*else2*)
                  key = ShiftFactor ;
		  If[ key =!= Infinity ,
		      (*then*)
                      CautionSolve[
			 IMultiplier[1/(wF+key),wF,var,x],var] (*caut*)
			 /.var:>y[x]
                      ,
		      (*else3*)
		      key =  ScaleFactor ;
		      If[ key =!= Infinity ,  (* scalable equation clause  *)
			  (*then*)            (*      may become obsolete  *)
			  CautionSolve[
			     IMultiplier[1/(var - key x wF),wF,var,x],var]
			     /.var:>y[x]
                          ,
			  (*else4*)
                          key = DilatationFactor;
					      (* dilatation test passed    *)
			  If[ key =!= Infinity,
			      (*then*)
                              CautionSolve[
				 IMultiplier[1/(var wF + key x),wF,var,x],
				 var ]
				 /.var:>y[x]
                              ,
			      (*else5*)
			      key = ScaleFactors;
			      If[ key =!= Infinity,
				  (*then*)
				  CautionSolve[
				  IMultiplier[1/(key[[2]] var + key[[3]] -
				     ( x + key[[1]]) wF),wF,var,x],
				     var ]
				     /.var:>y[x]
                                  ,
				  (*else6*)
			          key = DeepThought[ wF , var , x ] ;
			          If[ (key =!= Infinity) (*&&*) ,
				      (*then*)
				      key/.var:>y[x],
				      (*else*)
				      key = DeeperThought;
				      If[ key =!= Infinity,
				          (*then*)
				          key,
				          (*else*)
                                          key = DeepThought[1/wF, x , var ] ;
				          If[ key =!= Infinity,
					      (*then*)
				              (key/.var:>y)/.x:>x[y],
				              (*else*)
				              DSFail
					      (* the ultimate failure *)
                                            ] (*fi9*)
                                        ] (*fi8*)
                                    ] (*fi7*)
                                ] (*fi6*)
                            ] (*fi5*)
                        ] (*fi4*)
                    ] (*fi3*)
                ] (*fi2*) (*this was the if starting with Sep Var test    *)
           ]  (*fi1*)(*this was the if starting with the "no var" condition*)

      ] (* Module ends here *)
;
(****************** further front end considerations ***********************)

DS[ G_. y_'[x_] + H_. == 0/;(FreeQ[G,y'[x]] && FreeQ[H,y'[x]]) , y_[x_] ] :=
    Module[ {var, w, wG, wH, K, g },
	    CautionSolve [ eq_, term_ ] :=
	       Module [ {ws} , ws = Solve[ eq/.
		  { c_. Log[a_]+d_. Log[b_] :> Log [ a^c b^d ] }, term ] ;
		  If [ MatchQ[ws,{}],SolveFail[ eq, term ] , ws ] ]
		  (*eludom*) ;
	    wG = G/.y[x]:>var;
	    wH = H/.y[x]:>var;
	    K  = Simplify[Expand[D[wG,x]-D[wH,var]]];
	    If[ K  === 0,
		(*then*)
		w = Integrate[wG,var];
		CautionSolve[ w + Integrate[wH-D[w,x],x]==C[1],var]/.var:>y[x]
		,
		(*else*)
		g = Simplify[ K/wG];
		If[ FreeQ[g,var],
		    DS[Collect[Expand[Exp[-Integrate[g,x]](G y'[x] + H)],y'[x]]
		       ==0,y[x]]
		    ,
		    (*else*)
		    g = Simplify[ K/wH];
		    If[ FreeQ[g,x],
			DS[Collect[Expand[Exp[Integrate[g,var]/.var:>y[x]] *
			   (G y'[x] + H)],y'[x]]==0,y[x]]
                        ,
			(*else*)
		        DS[y'[x]==-H/G,y[x]]
                      ] (*fi*)
                   ] (*fi*)
               ]
           ] (*eludom*)

;
ClairotQ[ F_,y_[x_] ] :=
    Module [ { var , var1 , wF },
	     wF=F/.{y[x]:>var,y'[x]:>var1};
	     If[ FreeQ[wF,x] && FreeQ[wF,var],
		 (*then*)
		 False,
		 (*else*)
                 Simplify[Expand[
                 D[wF,x] + var1 D[wF,var]
		     ]] === 0
               ] (*fi*)
           ] (*eludom*)
;
DS[ F_==0 /; ClairotQ[F,y[x]] , y_[x_]] :=
    Module[
	   {u, wF, wS},
	   wF = F/.{y[x]:>u+x y'[x] } ;
	   wF = Simplify[Expand[wF]];
	   (* Must depend only on y'[x] and u now *)
	   wF = wF/.{y'[x]:>C[1]} ;
	   wS = Solve[wF==0,u];
	   If [ MatchQ[wS,{__}],
		wS/.{(u->r_):>(y[x]->C[1] x + r)},
		(*else*)
		DSFail
              ] (*fi*)
          ] (*eludom*)
;

(* Update as of 2/15/92: provide for equations not resolved *)
(* w.r.t. the first derivative                              *)

DS[lhs_==rhs_ ,y_[x_]] :=

    Module[ {w,DSList} ,
	    DSList[{}] := {};
	    DSList[{{y'[x]->term_},L___}] :=
	       Union[DSolve[y'[x]==term,y[x],x],DSList[{L}]];
            DSList[other_] := DSFail;
	    w=Simplify[Expand[lhs-rhs]];
	    If[ FreeQ[w,y'[x]],
		Solve[w==0,y[x]],
		If[ D[D[w,y'[x]],y'[x]] === 0 ,
		    DS[Collect[w,y'[x]] == 0 , y[x] ],
		    (*else*)
		    DSList[Solve[w==0,y'[x]]]
                  ]
              ]
          ]
;
DS[pro_ , y_[x_], x_ ] := DS[pro,y[x]]
;
DS[{lhs_==rhs_,c_. y_[a_]==b_},y_[x_]]
	    /; (FreeQ[{a,b,c},x] && FreeQ[{a,b,c},y] ) :=
            Module[ { Instantiate , memo },
	    Instantiate[{y[x]->term_},aa_,bb_] :=
	       Module[ { s },
		       s = Solve[bb ==(term/.x:>aa),C[1]];
		       If[ MatchQ[s,{}],
			   (*then*)
			   Message[DSolve::dsing, y[a]==b, a];
			    {},
                           (*else*)
			   If[ MatchQ[ s , {L__} ],
			       (*then*)
			       {y[x]->term}/.s
			       ,
			       (*else*)
			       {{y[x]->term},s}
                             ] (*fi*)
                         ] (*fi*)
                     ] (*eludom*)
		     ;
             Off[Solve::eqf];

	     Instantiate[Solve[other_,y[x]],aa_,bb_] :=
		Module[ { s }.
			s = Solve[(other/.{y[x]:>bb,x:>aa}),C[1]];
			If[ MatchQ[ s , {L__} ],
			    Solve[other,y[x]]/.s
			    ,
			    {Solve[other,y[x]],s}
                          ]
                     ]
		     ;
		     On[Solve::eqf];

	     Instantiate[{},aa_,bb_] := {};

	     Instantiate[{first_,rest___},aa_,bb_] :=
		 Union[Instantiate[first,aa,bb],Instantiate[{rest},aa,bb]];

	     Instantiate[other_,aa_,bb_]:= other ;

              memo = DS[lhs==rhs,y[x]];
	     Instantiate[memo,a,b/c]
          ] (*eludom*)
	  ;

DS[{c_. y_[a_]==b_,lhs_==rhs_},y_[x_]]
	    /; (FreeQ[{a,b,c},x] && FreeQ[{a,b,c},y] ) :=
	    DS[{lhs==rhs,c y[a] == b},y[x]]
;
DS[ pro_ , y_ , x_ ] /; (MatchQ[pro,lhs_==rhs_] ||
			MatchQ[pro,{lhs_==rhs_,c_. y[a_] == b_ }] ||
			MatchQ[pro,{c_. y[a_] == b_ , lhs_==rhs_}] ) :=

    Module[ {Purify, memo},
	    Purify[{u_[z_]->term_}] :=
		   {u->(Function[term]/.z:>#)};
            Purify[S_] :=  (S/.y[x]:>y)/.x:># ;
            memo=DS[pro,y[x]];
	    If[ MatchQ[memo, {___}],
		Map[Purify,memo],
		(*else*)
		Purify[memo]
              ] (*fi*)
	  ] (*eludom*)
          ;

DS[otherwise___]/;True:=DSFail;

(*********************** The main interface part ******************)
   memo=DS[prob,uu,t];
   memo /; FreeQ[memo, DSFail]
]/; Applicable1[prob, uu, t]

Applicable1[GG_==HH_,uu_,t_] := Applicable11[GG,HH,uu,t]

Applicable1[{GG_==HH_,cc_. yy_[aa_]==bb_},uu_,t_] := Applicable11[GG,HH,uu,t]

Applicable1[{cc_. yy_[aa_]==bb_,GG_==HH_},uu_,t_] := Applicable11[GG,HH,uu,t]

Applicable1[other_]:= False

Applicable11[GG_, HH_, uu_, t_]:= Module[ {h},
	      If [ MatchQ[uu , y_[t] ] ,
		   (*then*)
		   h=Head[uu],
		   (*else*)
		   h=uu
                 ] (*fi*)
		 ;
	      If [ Union[Cases[GG==HH,Derivative[_][h][t],-1 ]] =!= {h'[t]},
                   (* then it is not a first- order equation *)
		    False,
		   (*else*)
		   ( D[D[GG-HH,h[t]],h[t]] =!= 0 ||
		       D[D[GG-HH,h'[t]],h'[t]] =!= 0 ||
		       D[D[GG-HH,h[t]],h'[t]] =!= 0 )
                 ] (*fi*)
           ]

Unprotect[D]

D[integrate[f_,{var_,lower_,upper_}],x_]:= integrate[D[f,x],{var,lower,upper}] +
(f/.{var:>upper}) D[upper,x]

D[integrate[f_,u_],x_] :=
       Module[ {var},
       If [ FreeQ[u,x],
	    Integrate[ D[f,x], u] ,
	    (* else *)
            Integrate[ D[f/.{u:>var},x],{var,0,u}] +
	    f D[u,x]
          ] (*fi*)
	  ] (*eludom*)

D[Integrate[f_,{var_,lower_,upper_}],x_]:= Integrate[D[f,x],{var,lower,upper}] +
					  (f/.{var:>upper}) D[upper,x]

Protect[D]


DSolve[prob_,uu_,t_] :=
   Module[ {DS2, memo, DSNormalize, DSFail },

DS2[y_''[x_]==F_,y_[x_]] :=

     Module [ { CautionSub , CautionIntegrate , CrazySolve ,
                IMultiplier , ScaleSpecial , ScaleSpecial1 ,
                SepVarSpecial , ShiftFactor ,
                Scale2Special , ClairotwoQ , EasyTotal , FirstObstacleQ ,
                SecondObstacleQ , SepVar , key , wF, var , var1 , u , z ,
                memo , LSFQ , LSFIQ , LS2F , LSV , wwF ,
                RemoveFirstObstacle , SeparateVariables } ,

(*************** Here goes the 2nd argument of this Module ***************)

   CautionSub[expr_,lhs_,rhs_]:=
    Module[ { w , CS },
	    CS[f_,u_,l_,r_] :=
		  If [ u =!= (u/.l:>r),
		       (*then*)
		       Integrate[ (f/.u:>var) , {var,0,u} ],
		       (*else*)
		       Integrate[ f , u ]
                     ] ;

	    w = expr/.Integrate[f_,u_]:> CS[f,u,lhs,rhs];
            (w/.{lhs:>rhs})
          ]  (*eludom*) ;


   CautionIntegrate[ expr_ , yy_ ] :=
     Module [ {v, memo },
	      memo = Integrate[expr,yy];
	      If[ FreeQ[ memo , Integrate[f_,u_] ],
		   (*then*)
		   memo,
		   (*else*)
		   Integrate[(expr/.{yy:>v}),{v,0,yy}]
                ] (*fi*)
            ] (*eludom*) ;


   CrazySolve [ eq_, term_ ] :=
    Module [ {ws} ,
              If[ FreeQ[eq,Integrate[f_,u_]] && FreeQ[eq,integrate[f_,u_]],
	      Off[Solve::dinv];
	      Off[General::stop];
	      ws = Solve[ eq
	      , term ] ;
	      On[Solve::dinv];
	      On[General::stop];
	      If [ MatchQ[ws,{}],
	      SolveFail[ eq, term ] , ws ]
              ,
              {eq}
	      ] ] ;

   IMultiplier [ m_ , FF_ , yy_[xx_] ] :=
       Module [ { w } ,
		w= CautionIntegrate [ m,yy[xx]] ;
		CrazySolve [
		     (( w - Integrate [ Simplify[
		     D[ w/.yy[xx]:>var , xx ] +
		     (FF m )/.yy[xx]:>var ] , xx ] )
		     /.var:>yy[xx]) == C[1] , yy[xx]
		     ]
              ]  ;

    ScaleSpecial :=
	 Module[ { Fx, Fy, Fxx, Fxy, Fyy, Fxyy, Fyyy, a, Den  } ,
		 Fx = D[wF,x];
		 Fy = D[wF,var1];
		 Fxx = D[Fx,x];
		 Fxy = D[Fx,var1];
		 Fyy = D[Fy,var1];
		 Fxyy = D[Fyy,x];
		 Fyyy = D[Fyy,var1];
		 Den=Simplify[
    (2*Fxyy*wF - 2*Fxyy*Fy*var1 + Fx*Fyyy*var1 + Fxyy*Fyy*var1^2 - Fxy*Fyyy*var1^2)
			     ];
             If[ Den===0,
		 Infinity,
		 a= Simplify[
    (-2*Fyy*wF + 2*Fy*Fyy*var1 - Fyyy*wF*var1 - Fyy^2*var1^2 + Fy*Fyyy*var1^2)
		 /Den
		];
		If [ FreeQ[a,var1],
		     If [ Simplify[Expand[
			     (D[a,x]-1) var1 Fy - a Fx + (1-2 D[a,x] ) wF-
			     D[D[a,x],x] var1 ]]===0,
                          a,
			  Infinity
                        ] 
			,
                     Infinity
                   ] 
               ] 
               ]  ;


   ScaleSpecial1 :=
   Module[ { Fy, Fy1, Fyy1, Fy1y1, Fyy1y1, Fy1y1y1, Den, b },
	   Fy=D[wF,var];
	   Fy1=D[wF,var1];
	   Fyy1=D[Fy,var1];
	   Fy1y1=D[Fy1,var1];
	   Fyy1y1=D[Fyy1,var1];
	   Fy1y1y1=D[Fy1y1,var1];
           Den = Simplify[
  (2*Fyy1*wF - 2*Fyy1*Fy1*var1 - 2*Fyy1y1*wF*var1 + 2*Fyy1y1*Fy1*var1^2 +
   Fyy1*Fy1y1*var1^2 - 2*Fy*Fy1y1y1*var1^2 - Fyy1y1*Fy1y1*var1^3 +
   Fyy1*Fy1y1y1*var1^3)
			 ] ;
           If [ Den === 0,
		Infinity,        
                b = Simplify[
 (-2*Fy1*wF + 2*Fy1^2*var1 + 2*Fy1y1*wF*var1 - 3*Fy1*Fy1y1*var1^2 +
     2*Fy1y1y1*wF*var1^2 + Fy1y1^2*var1^3 - Fy1*Fy1y1y1*var1^3)/
			     Den ];
                If [ FreeQ[b,var1],
		     If [ Simplify[Expand[
			     (D[b,var]-1) var1 Fy1 + b Fy + (2 - D[b,var]) *
			     wF - D[D[b,var],var] var1^2 ]] === 0,
			  b,
			  Infinity
                        ]
			,
                     Infinity
                  ] 
             ] 
        ]  ;

  SepVarSpecial :=

     Module [ { w, Den } ,
              Den=Simplify[Expand[var1^2 D[D[wF,var1],var1] -
				  2 var1 D[wF,var1] + 2 wF ] ] ;
              If [ Den === 0,
		   Infinity,      
                   w = Simplify[Expand[
		       2 D[wF,var] - var1 D[D[wF,var],var1] ]/Den] ;
                   If [ FreeQ[w,var1],
			If[ Simplify[Expand[
			       D[wF,var] + w var1 D[wF,var1] -
			       w wF - ( D[w,var] + w^2 ) var1^2 ]] === 0,
                            Exp[Integrate[w,var]],
			     Infinity
                          ] 
			  ,
			Infinity
                      ] 
                 ] 
            ]  ;


  ShiftFactor  :=
     Module [ { w } ,
	      w = Simplify[ D[wF,x] / D[wF,var] ] ;
	      If [ FreeQ[w,var] && FreeQ[w,x] && FreeQ[w,var1] ,
		    w,
		    Infinity ]
            ] ;



   Scale2Factor :=
    Module [ { b },

	    b = Simplify[Expand[
		( 2 wF + x D[wF,x] - var1 D[wF,var1])/
		( wF - var D[wF,var] - var1 D[wF,var1] )
		]];

           If [ FreeQ[b,x] && FreeQ[b,var] && FreeQ[b,var1],
		b,
		Infinity
              ] 
           ]  ;

  ClairotwoQ :=
    Module[ { },
            Simplify[Expand[D[wF,x]+var1 D[wF,var]]] === 0
          ]  ;


   EasyTotal :=
    Module[ {l, Fx, Fy, Fyy1, Fxy1, eq, Den},
	    Fx=D[wF,x];
	    Fy=D[wF,var];
	    Den=Simplify[Expand[(Fx+var1 Fy)]];
	      Fy1=D[wF,var1];
	      Fxy1=D[Fx,var1];
	      Fyy1=D[Fy,var1];
	    l= Simplify[Expand[(Fy - Fyy1 var1 - Fxy1)/Den]];
	    If [ FreeQ[l,x] && FreeQ[l,var] ,
		 If [ Simplify[Expand[
		      wF (D[l,var1] + l^2) + 2 Fy1 l + D[Fy1,var1]
		      ]] === 0,
		      l = Exp[Integrate[l,var1]];
		      l/.var1:>y'[x],
		      Infinity
                    ] 
		    ,
		 Infinity
               ] 
          ]  ;

   FirstObstacleQ :=

    Module[ {},
	    Simplify[D[D[wF,x],var]]===0
          ]  ;

   RemoveFirstObstacle :=
    Module[ {l, Fx, Fy, Fyy1, Fxy1, eq },
            Fx=D[wF,x];
	    Fy=D[wF,var];
	    Fy1=D[wF,var1];
	    Fxy1=D[Fx,var1];
	    Fyy1=D[Fy,var1];
            l= Simplify[Expand[(Fy - Fyy1 var1 - Fxy1)/(Fx+var1 Fy)]];
	    If [ FreeQ[l,x] && FreeQ[l,var] ,
		 If [ Simplify[Expand[
		      wF (D[l,var1] + l^2) + 2 Fy1 l + D[Fy1,var1]
		      ]] === 0,
		      l = Exp[Integrate[l,var1]];
                      eq= Integrate[(l/.var1:>y'[x]) (y''[x] - F ) , x]
                           == C[2];
		      DSolve[eq,y[x],x]
		      ,
                      DSFail
                    ] ,
		    DSFail
               ] 
           ] 
           ;

   SecondObstacleQ :=

    Module[ {},
	    Simplify[var1 D[D[wF,x],var1]-D[wF,x]]===0
          ]  ;


   SepVar :=

    Module[ { g1, Den, Num },
	    Den=Simplify[Expand[D[wF,x] - var1*D[D[wF,x],var1] ]];
	    Num = Simplify[Expand[D[D[wF,x],var]]];
	    If [ Den === 0 ,
		 If [ Num =!= 0 ,
		      Infinity,
		      Pending
                    ] 
                 ,
                 g1 = Simplify[Num/Den] ;
	         If [ FreeQ[g1,x] && FreeQ[g1,var1],
                      If[ Simplify[Expand[
		          D[wF,var]/g1 + var1 D[wF,var1] - wF -
		          (D[g1,var]/g1+g1) var1^2 ]] === 0,
		          Exp[Integrate[g1,var]]/.var:>y[x],
		          Infinity
                        ] 
		        ,
		        Infinity
                     ] 
                ] 

          ]  ;

   SeparateVariables :=
	   Module[ { },
		   wF = PowerExpand[Simplify[Expand[(F/.y'[x]:>
                                                    (var1[x] key ))/key -
				                    D[key,y[x]] var1[x]^2
                                                    ]]];
		   If [ FreeQ[wF,y[x]],
			memo = DSolve[var1'[x]==wF,var1[x],x];
			If [ MatchQ[Head[memo],DSolve],
			     DSFail,
                             LSV[memo]
                          ] 
			,
			DSFail
                      ] 
                ] 
                ;

          LSFQ[{}]:={};
	  LSFQ[{var1[x]->term_}]:= {y[x]->Integrate[term,x]+C[2]}; 
	  LSFQ[{first_,L___}] :=
                 Prepend[LSFQ[{L}],LSFQ[first]];
          LSFQ[other_]:={other/.{var1[x]:>y'[x]}} ;

          Off[Union::heads];
          LSFIQ[{}]:={};
	  LSFIQ[{var1[z]->term_}]:=
                 CrazySolve[(Integrate[1/term,z]/.z:>y[x])==x+C[2],
                       y[x]] ;
          LSFIQ[{first_,L__}]:= Union[LSFIQ[first],LSFIQ[{L}]];
          LSFIQ[{only_}] := LSFIQ[only] ;
          LSFIQ[other_]:={other/.{var1[z]:>y'[x], z:>y[x]}} ;

          LS2F[{u[z]->term_}] :=
		 IMultiplier[
		    1/Simplify[ key y[x] -
		    x^key CautionSub[term,z,y[x]/x^key]],
		    x^(key-1) CautionSub[term,z,y[x]/x^key], y[x] ];
          LS2F[{}] := {};
          LS2F[{first_,rest__}] :=
		 Union[LS2F[first],LS2F[{rest}]];
          LS2F[{only_}] := 
	  LS2F[only]
	  ;		  
          LS2F[other_] := 
	  {other/.{u[z]:>(y'[x] x^(1-key)),z:>(y[x]/x^key)}} 
	  ;

          LSV[ {var1[x]->term_} ] :=
				 CrazySolve[ Integrate[1/key,y[x]] ==
				    Integrate[term,x] + C[2] , y[x] ];
          LSV[{}] := {};
          LSV[{first_,rest__}] :=
                Union[LSV[first],LSV[{rest}]];
          LSV[{only_}]:=LSV[only] ;

          LSV[other_] := {other/.
				  {  var1[x] :> (y'[x]/key) }} ;
          On[Union::heads];
(************** the main body of the module starts here ****************)
          Off[Union::heads];
wF = (F/.y[x]:>var)/.y'[x]:>var1 ;  
If[ FreeQ[wF,var] ,
    key= ScaleSpecial;
    If [ key=!=Infinity,
	 memo= IMultiplier[1/((1-D[key,x]) var1[x] - key (wF/.{var1:>var1[x]})),
			   (wF/.{var1:>var1[x]}),
			   var1[x]]
         ,
         memo=DSolve[var1'[x]==(wF/.var1:>var1[x]),var1[x],x];
        ]  ;
     If[ MatchQ[Head[memo],DSolve]
         ,
         DSFail,
         LSFQ[memo]
        ] 
      ,
      (*1else*)
      If[ FreeQ[wF,x],
	  wwF=((F/.y[x]:>z)/.y'[x]:>var1[z])/var1[z];
	  key= ScaleSpecial1;
	  If [ key=!=Infinity,
               key=key/.var:>z;
               memo= IMultiplier[1/((D[key,z]-1) var1[z] -key wwF), wwF,
		                 var1[z]]
               ,
	       key=SepVarSpecial;
	       If [ key=!=Infinity,
		    key=key/.var:>z;
                    memo= IMultiplier[1/(D[key,z] var1[z] - key wwF), wwF,
		                      var1[z]],
		    memo=DSolve[var1'[z]==wwF,var1[z],z]
                  ] 
             ] 
	     ;
          If[ MatchQ[Head[memo],DSolve]
	      ,
	      DSFail,
              LSFIQ[memo]
             ]
          ,
          (*2else*) 
          key =  Simplify[Expand[ F x y[x]^2 ]] ;
          If [ FreeQ[ key , x ] ,
	       {Solve[Log[-1 + 2 C[1] y[x] / x + 2 Sqrt[C[1]]/x
			Sqrt[C[1] y[x]^2 - x y[x] ] ] /(2 C[1]^(3/2)) +
                        Sqrt[C[1] y[x]^2 - x y[x] ] / (C[1] x) ==
			(2 key)/ x + C[2] , y[x] ],
                     Solve[Log[-1 + 2 C[1] y[x] / x + 2 Sqrt[C[1]]/x
			Sqrt[C[1] y[x]^2 - x y[x] ] ] /(2 C[1]^(3/2)) +
			Sqrt[C[1] y[x]^2 - x y[x] ] / (C[1] x) ==
			- (2 key)/ x + C[2] , y[x] ] }
               ,
               (*3else*) 
               key:=TotalDerivative[F,y[x]] ;
               If[ key=!=Infinity,
         	   memo = DSolve[y'[x] == key + C[2], y[x],x];
	           If [ MatchQ[Head[memo],DSolve],
		        DSFail,
		        memo
                      ]
                      ,
                   (*4else*)
                   If[Simplify[F - D[F,y[x]] y[x] - D[F,y'[x]] y'[x] ] === 0
                      ,
	              key = ( F/y[x] )/.{y'[x]:>u[x] y[x]};
	              key = PowerExpand[Simplify[Expand[key]]];
	              If [ FreeQ[key,y[x]],
		           memo:=DSolve[u'[x]==key - u[x]^2, u[x],x];
		           If [ MatchQ[ Head[memo] , DSolve ],
		                DSFail,
                                memo/.
		                {{u[x]->term_}:>
                                 {y[x]->C[2] Exp[Integrate[term,x]]},
		                 u:>y'}
                              ] 
		          ,
		          DSFail
                         ] 
                      ,
                      (*5else*)
                      key = ShiftFactor;
                      If[ key =!=Infinity,
	                  wwF = Simplify[F/.{y[x]:>(z-w x),y'[x]:>var1[z]}] ;
	                  If [ FreeQ[wwF,x],
                               memo=DSolve[ (var1[z]+key) var1'[z] == (wwF) ,
                                    var1[z],z];
		               If[ MatchQ[Head[memo],DSolve],
		                   DSFail,
                                   memo=memo/.C[1]:>C[2];
		                   memo/.{{var1[z]->term_} :>
		                   IMultiplier[1/(key+(term/.z:>
                                               (y[x]+key x))),
				               (term/.z:>(y[x]+key x)),
                                               y[x] ] ,
                                               z :> (y[x]+key x) , var1:>y'}
                                 ] 
		               ,
		               DSFail
                             ] 
                          ,
                          (*6else*)
                          key = Scale2Factor;
                          If[ key=!=Infinity,
	                      wwF = PowerExpand[Simplify[Expand[
		              ((x^(2-key)F+(1-key) u[z])/(u[z]-key z))/.
		              {y[x]:>(z x^key),y'[x]:>(u[z] x^(key-1) )}
                                		  ]]];
                              If [ FreeQ[wwF,x] ,
		                   memo = DSolve[ u'[z] == wwF , u[z],z ];
		                   If [ MatchQ [ Head[memo], DSolve ]
			                ,
			                DSFail,
                                        memo=memo/.C[1]:>C[2];
			                LS2F[memo]
                                      ] 
		                   ,
		                   DSFail
                                 ] 
                              ,
                              (*7else*)
                              If[ ClairotwoQ,
	                          wwF=F/.{y'[x]:>z,y[x]:>(u[z]- z u'[z]),
                                          x:>(-u'[z])};
	                          If [ FreeQ[wwF,x],
		                       memo := DS2[u''[z]==-1/wwF,u[z]]
                          	       ;
		                       memo/.
		                       {{u[z]->term_}:>
                                        (ParametricForm[x->-D[term,z],
						y-> term - z D[term,z]]
						/.z:>DSolve`z),
                                        Solve[aeq_,u[z]]:>
                                        (ParametricForm[x ->
                                        -((u'[z]/.Solve[D[aeq,z],u'[z]])[[1]]),
                                        Solve[(aeq/.{u[z]:>(y - x z)}),y]]/.
                                        z:>Dsolve`z)
                                       } 
		                       ,
		                       DSFail
                                     ] 
                                  ,
                                  (*8else*)
                                  key =EasyTotal;
                                  If[ key =!=Infinity,
	                              DSolve[Integrate[key y''[x]-key F,x]==
                                         C[2],y[x],x]
                                      ,
                                      (*9else*)
                                      If [ FirstObstacleQ,
                                           RemoveFirstObstacle,
                                           (*10else*)
                                           If[ SecondObstacleQ,
                                               DSFail,
                                               (*11else*)
                                               key = SepVar;
                                               If[ key =!=Infinity,
                                                   SeparateVariables,
                                                   DSFail
                                                 ] (*fi11*)
                                             ] (*fi10*)
                                         ] (*fi9*)
                                    ] (*fi8*)
                                ] (*fi7*)
                            ] (*fi6*)
                        ] (*fi5*)
                     ] (*fi4*)
                 ] (*fi3*)
             ] (*fi2*)
        ] (*fi1*)
  ] (*fi*)
] (*eludom*)  (* The main , hardworking DS2 module ends here *)
;
TotalDerivative[F_,y_[x_]]:=
    Module[ { IntSuccess , memo, Casualties },
	    IntSuccess[{}] := True;
	    IntSuccess[{Integrate[f_,tt_], L___}] :=
	       FreeQ[f,y'[x]] && FreeQ[f,y[x]] && IntSuccess[{L}] ;
            memo= Integrate[F,x];
	    Casualties = If[ MatchQ[memo,Integrate[g_,z_]],
			     {memo},
			     (*else*)
			     Cases[memo,Integrate[g_,z_],Infinity]
                           ] (*fi*)
                           ;
	    If [ IntSuccess[Casualties],
		 memo,
		 (*else*)
		 Infinity
               ] (*fi*)
          ] (*eludom*)
;

DS2[(G_==F_) ,y_[x_]] :=
     Module[ {memo },
	     memo = TotalDerivative[G-F,y[x]];
	     memo = DSolve[ memo == C[2], y[x],x];
	     If [ MatchQ[Head[memo],DSolve],
		  DSFail,
		  (*else*)
		  memo
                ]
           ] (*eludom*)  /; TotalDerivative[G-F,y[x]]=!=Infinity
;
DS2[lhs_==rhs_ ,y_[x_]] :=
    Module[ {w,DSList} ,
	    DSList[{}]:= {};
	    DSList[{{y''[x]->first_},rest___}] :=
	       Union[DSolve[y''[x]==first,y[x],x],DSList[{rest}]];
            DSList[other_]:=DSFail;
	    w=Simplify[Expand[lhs-rhs]];
	    If[ FreeQ[w,y''[x]],
		DSolve[w==0,y[x],x],
		(*else*)
		DSList[Solve[w==0,y''[x]]]
              ] (*fi*)
          ] (*eludom*)
;

DS2[pro_ , y_[x_], x_ ] := DS2[pro,y[x]]
;

DS2[{lhs_==rhs_,c0_. y_[a_] == b0_ , c1_. y_'[a_] == b1_} , y_[x_] ]
    /; (FreeQ[{a,b0,b1,c0,c1},x] && FreeQ[{a,b0,b1,c0,c1},y] ) :=
    Module[ { Instantiate , memo },
	    Instantiate[{y[x]->term_},aa_,bb0_,bb1_] :=
	       Module[ { s },
		       s = Solve[{bb0==(term/.x:>aa),
				  bb1==(D[term,x]/.x:>aa)},
				  {C[1],C[2]}
                                ];
                       If[ MatchQ[s,{}],
			   (*then*)
			   Message[DSolve::dsing,
                               {c0 y[a] == b0, c1 y'[a] == b1}, a];
			   {},
			   (*else*)
			   If[ MatchQ[ s , {L__} ],
			       (*then*)
			       {y[x]->term}/.s
			       ,
			       (*else*)
			       {{y[x]->term},s}
                             ] (*fi*)
                         ] (*fi*)
                     ] (*eludom*)
		     ;

	     Instantiate[{},aa_,bb0_,bb1_] := {};

	     Instantiate[{first_,rest___},aa_,bb0_,bb1_] :=
		 Union[Instantiate[first,aa,bb0,bb1],
			Instantiate[{rest},aa,bb0,bb1]];

             Instantiate[other_,aa_,bb0_,bb1_]:= {other, c0 y[a] == b0,
						   c1 y'[a] == b1 };

	     memo = DS2[lhs==rhs,y[x]];
	     Instantiate[memo,a,b0/c0,b1/c1]
          ] (*eludom*)
	  ;

DS2[{lhs_==rhs_,c1_. y_'[a_] == b1_,c0_. y_[a_] == b0_}, y_[x_] ]
    /; (FreeQ[{a,b0,b1,c0,c1},x] && FreeQ[{a,b0,b1,c0,c1},y] ) :=
    DS2[{lhs==rhs,c0 y[a] == b0,c1 y'[a] == b1},y[x]] ;

DS2[{c0_. y_[a_] == b0_,lhs_==rhs_,c1_. y_'[a_] == b1_}, y_[x_] ]
    /; (FreeQ[{a,b0,b1,c0,c1},x] && FreeQ[{a,b0,b1,c0,c1},y] ) :=
    DS2[{lhs==rhs,c0 y[a] == b0,c1 y'[a] == b1},y[x]] ;

DS2[{c0_. y_[a_] == b0_,c1_. y_'[a_] == b1_,lhs_==rhs_}, y_[x_] ]
    /; (FreeQ[{a,b0,b1,c0,c1},x] && FreeQ[{a,b0,b1,c0,c1},y] ) :=
    DS2[{lhs==rhs,c0 y[a] == b0,c1 y'[a] == b1},y[x]] ;

DS2[{c1_. y_'[a_] == b1_,lhs_==rhs_,c0_. y_[a_] == b0_}, y_[x_] ]
    /; (FreeQ[{a,b0,b1,c0,c1},x] && FreeQ[{a,b0,b1,c0,c1},y] ) :=
    DS2[{lhs==rhs,c0 y[a] == b0,c1 y'[a] == b1},y[x]] ;

DS2[{c1_. y_'[a_] == b1_,c0_. y_[a_] == b0_,lhs_==rhs_}, y_[x_] ]
    /; (FreeQ[{a,b0,b1,c0,c1},x] && FreeQ[{a,b0,b1,c0,c1},y] ) :=
    DS2[{lhs==rhs,c0 y[a] == b0,c1 y'[a] == b1},y[x]] ;


DS2[pro_ , y_ , x_ ] 
     /; (MatchQ[pro,lhs_==rhs_] ||
	 MatchQ[pro,{lhs_==rhs_, c0_. y[a_] == b0_ , c1_. y'[a_] == b1_ }] ||
	 MatchQ[pro,{lhs_==rhs_, c1_. y'[a_] == b1_ , c0_. y[a_] == b0_ }] ||
	  MatchQ[pro,{c0_. y[a_] == b0_ ,lhs_==rhs_, c1_. y'[a_] == b1_ }] ||
	  MatchQ[pro,{c0_. y[a_] == b0_ ,c1_. y'[a_] == b1_,lhs_==rhs_ }] ||
	  MatchQ[pro,{c1_. y'[a_] == b1_,lhs_==rhs_, c0_. y[a_] == b0_ }] ||
	  MatchQ[pro,{c1_. y'[a_] == b1_,c0_. y[a_] == b0_ ,lhs_==rhs_ }] )
:=
    Module[ {Purify, memo},
	    Purify[{u_[z_]->term_}] :=
		   {u->(Function[term]/.z:>#)};
            Purify[S_] := (S/.y[x]:>y)/.x:># ;
	    memo=DSNormalize[DS2[pro,y[x]]];
            If[ MatchQ[memo, {___}],
		Map[Purify,memo],
		(*else*)
		Purify[memo]
              ] (*fi*)
          ] (*eludom*)
;
DS2[otherwise___]/;True := DSFail
;
DSNormalize[{}] := {}
;
DSNormalize[{y_[x_]->term_,rest___}] :=
     Prepend[ DSNormalize[{rest}],{y[x]->term}]
;
DSNormalize[{y_->func_,rest___}] :=
     Prepend[ DSNormalize[{rest}],{y->func}]
;
DSNormalize[Union[{},S_]]:= S 
;
DSNormalize[other_] := other
;
(********************* The main interface part ***********************)
   memo= DSNormalize[DS2[prob,uu,t]];
   memo /; FreeQ[memo, DSFail]
] /;Applicable2[prob, uu, t]

Applicable2[GG_==HH_,uu_,t_] := Applicable22[GG,HH,uu,t]

Applicable2[{GG_==HH_,c0_. y_[a_] == b0_ , c1_. y_'[a_] == b1_ },uu_,t_] :=
	   Applicable22[GG,HH,uu,t]

Applicable2[{GG_==HH_,c1_. y_'[a_] == b1_ ,c0_. y_[a_] == b0_ },uu_,t_] :=
	   Applicable22[GG,HH,uu,t]

Applicable2[{c0_. y_[a_] == b0_ ,GG_==HH_,c1_. y_'[a_] == b1_ },uu_,t_] :=
	   Applicable22[GG,HH,uu,t]

Applicable2[{c0_. y_[a_] == b0_ ,c1_. y_'[a_] == b1_ ,GG_==HH_},uu_,t_] :=
	   Applicable22[GG,HH,uu,t]

Applicable2[{c1_. y_'[a_] == b1_,GG_==HH_,c0_. y_[a_] == b0_},uu_,t_] :=
	   Applicable22[GG,HH,uu,t]

Applicable2[{c1_. y_'[a_] == b1_,c0_. y_[a_] == b0_,GG_==HH_},uu_,t_] :=
	   Applicable22[GG,HH,uu,t]

Applicable2[other_]:=False

Applicable22[GG_,HH_,uu_,t_] :=
   Module[{h,DerList},
	  If[ MatchQ[uu , y_[t] ] ,
	       (*then*)
	       h=Head[uu],
	       (*else*)
	       h=uu
            ] (*fi*)
	    ;
          DerList=Union[Cases[GG==HH,Derivative[_][h][t],-1 ]];
	  If[ (DerList=!={h'[t],h''[t]}) && (DerList=!={h''[t]}),
	      (* then it is not a second-order equation *)
	      False,
	      (*else*)
	      DerList=GG-HH;
	      D[D[DerList,h'[t]],h'[t]] =!= 0  ||
	      D[D[DerList,h[t]],h[t]] =!= 0  ||
	      D[D[DerList,h''[t]],h''[t]] =!= 0  ||
	      D[D[DerList,h'[t]],h[t]] =!= 0  ||
	      D[D[DerList,h'[t]],h''[t]] =!= 0  ||
	      D[D[DerList,h[t]],h''[t]] =!= 0
            ] (*fi*)
         ] (*eludom*)




 DSolve[ Equal[a_, b_], f_, x_ ] := 
   Module[ { answer },
	VariableC = olderdeg = True;
	countdsolve = 0; count=-1;
	answer = DSolveVictor[a==b,If[Head[f]===Symbol,f@@{x},f], x, 	
				DSolveConstants -> C];
	If[ FreeQ[answer, FailODE],
	    If[ Head[f]===Symbol,
		answer = {{f -> Function@@{ answer[[1,1,2]]/.x->#  } }}
	]];
	answer /; FreeQ[answer, FailODE]
   ]
   
 DSolve[ Equal[a_, b_], f_, x_, DSolveConstants -> g_ ] := 
   Module[ { answer },
	VariableC = olderdeg = True;
	countdsolve = 0; count=-1;
	answer = DSolveVictor[a==b,If[Head[f]===Symbol,f@@{x},f], x, 	
				DSolveConstants -> g];
	If[ FreeQ[answer, FailODE],
	    If[ Head[f]===Symbol,
		answer = {{f -> Function@@{ answer[[1,1,2]]/.x->#  } }}
	]];
	answer /; FreeQ[answer, FailODE]
   ]

 DSolve[ {Equal[a_, b_], cond___}, f_, x_ ] := 
   Module[ { answer, pure, p = If[Head[f]===Symbol,f,Head[f]] },
	VariableC = olderdeg = True;
	countdsolve = 0; count=-1;
	answer = DSolveVictor[a==b,If[Head[f]===Symbol,f@@{x},f], x, 	
				DSolveConstants -> C];
	If[ FreeQ[answer, FailODE],
	    pure = 
	    If[ Head[f]===Symbol,
		{{f -> Function@@{ answer[[1,1,2]]/.x->#} }},
		{{Head[f] -> Function@@{ answer[[1,1,2]]/.x->#} }} ];
	    p = Max[Cases[Expand[a-b], y_/;!FreeQ[y,Derivative[_][p][x]]]/.
	    	     w_. Derivative[n_][p][x] :> n];  
	    pure = Solve[ First[{cond}//.pure], Table[C[i], {i,p}] ];
	    answer = If[ Head[f]===Symbol,
			 {{f -> Function@@{
			 	   First[answer[[1,1,2]]/.x->#//.pure]} }},
			 {{f -> First[answer[[1,1,2]]//.pure]}}];				   
	];
	answer /; FreeQ[answer, FailODE]
   ]

 DSolve[ {Equal[a_, b_], cond___}, f_, x_, DSolveConstants -> g_  ] := 
   Module[ { answer, pure, p = If[Head[f]===Symbol,f,Head[f]] },
	VariableC = olderdeg = True;
	countdsolve = 0; count=-1;
	answer = DSolveVictor[a==b,If[Head[f]===Symbol,f@@{x},f], x, 	
				DSolveConstants -> g];
	If[ FreeQ[answer, FailODE],
	    pure = 
	    If[ Head[f]===Symbol,
		{{f -> Function@@{ answer[[1,1,2]]/.x->#} }},
		{{Head[f] -> Function@@{ answer[[1,1,2]]/.x->#} }} ];
	    p = Max[Cases[Expand[a-b], y_/;!FreeQ[y,Derivative[_][p][x]]]/.
	    	     w_. Derivative[n_][p][x] :> n];  
	    pure = Solve[ First[{cond}//.pure], Table[C[i], {i,p}] ];
	    answer = If[ Head[f]===Symbol,
			 {{f -> Function@@{
			 	   First[answer[[1,1,2]]/.x->#//.pure]} }},
			 {{f -> First[answer[[1,1,2]]//.pure]}}];				   
	];
	answer /; FreeQ[answer, FailODE]
   ]

(*****************************************************************************)
   
 DSolveVictor[ Equal[a_, b_], f_[x_], x_, DSolveConstants -> g_ ] := 
   Module[ { answer = Together[Expand[a - b]], p },
	Variable = Variable1 = True;
	countdsolve += 1;
	answer = 
	If[ !FreeQ[Denominator[ answer ], f[x]],
	    FailODE,
	    If[ MatchQ[Denominator[answer], Exp[s_/;!FreeQ[s,x]]],
		answer,
		Numerator[answer]
	   ]
	];
	answer = 
	If[ FreeQ[answer,FailODE],
	    p = Cases[answer, w_. Derivative[n_][_][x] ];
	    If[ !FreeQ[ Expand[answer - Plus@@p], Derivative ],
		FailODE,
		p = p/.w_. Derivative[n_][f][x] :> n;
		ODE[ answer, f[x], x, Max[p],
		     If[ FreeQ[answer, f[x]],
			 Min[p], 0], 
		     C ]
	    ],
	    FailODE
	]; 
	If[ count==-1, count = countdsolve-1, count-=1 ]; 
	If[ And@@(!FreeQ[answer, #]&/@{FailODE, DirectedInfinity[], 	
					Indeterminate}),
	    If[ countdsolve > 0 || i==0,
		countdsolve = 0;
		answer = 
		Block[ {Integrate,t},
		   Off[DSolve::dnim];
		   t=DSolve[ a==b,f[x],x,DSolveConstants -> g ];
		   On[DSolve::dnim];
		   If[ count==0, dsolve = True];
		   Block[ {DSolve},
			  If[ FreeQ[t,DSolve] ,
			      {{f[x]->ConstDelDSolve[t[[1,1,2]],x,Max[p],g]}}
			      ,
			      FailODE ]
		   ]
		],
		answer = FailODE
	   ]
	   ,
	   answer = If[ Head[answer]===List,
	   		{{f[x]->ConstDelDSolve[answer[[1,1,2]],x, Max[p], g]}},
			answer]
	];
	answer 
   ]
  
 DSolveVictor[ w__ ] := FailODE

 ODE[ __, max_?(Or[!IntegerQ[#],#<2]&), min_, w_ ] := 
	FailODE /; True

 ODE[ expr_, f_[x_], x_, max_, n_?(Or[!IntegerQ[#],#>0]&), w_ ] := 
   Module[ { newexpr, newf },
	newexpr = expr/.Derivative[k_][f][x] :> Derivative[k-n][newf][x];
	newexpr = DSolve[newexpr==0,newf[x],x,DSolveConstants -> w];
	newexpr = 
	Block[ {DSolve},
		If[ FreeQ[newexpr,DSolve],
		    newexpr[[1,1,2]],
		    FailODE ]
	];
	If[ FreeQ[newexpr, FailODE],
	    {{ f[x] ->
	    Nest[ Integrate[#,x]&, newexpr, n ] +
	    Plus@@Table[w[i] x^(i-max+n-1),{i,max-n+1,max}]
	    }},
	    FailODE
	]
   ] 

 ODE[ expr_, f_[x_], x_, 2, 0, w_ ] := 
   ODEsecond[Collect[Collect[Collect[ expr,f''[x]],f'[x]],f[x]],f[x],x,w];

 ODE[ expr_, f_[x_], x_, n_/;n>2, 0, w_ ] := ODEany[expr,f[x],x,n,w]

 ODE[ __ ] := FailODE

 ODEsecond[ expr_, f_[x_], x_, w_ ] :=
   Module[ { r },
     Off[ Power::infy, Infinity::indet ];
     r = Kummer[ expr, f[x], x, w];
     On[ Power::infy, Infinity::indet ];
     {{f[x] -> ConstDelDSolve[r[[1,1,2]],x,2,w] }} /;
     FreeQ[r,FailODE]
   ]

 ODEsecond[ a_. Derivative[2][f_][x_] + b_. Derivative[1][f_][x_] + c_,
	f_[x_], x_, w_ ] := 
   Module[ { newfun,answer },
	answer = (a f''[x] + b f'[x] + c)/.
	{
	Derivative[p_][f][x] :> D[E^(-b/(2a) x) newfun[x] ,{x,p} ],
	f[x] :> E^(-b/(2 a) x) newfun[x] 
	};
	answer = DSolveVictor[ E^(b/(2a) x) answer == 0, newfun[x], x, 
		 DSolveConstants->w ];
	answer = 
	If[ FreeQ[answer, FailODE],
	    {{ Rule[f[x], SSimplify[answer[[1,1,2]] E^(-b/(2a) x)]] }},
	    FailODE
	];
	answer /; FreeQ[answer, FailODE]
   ] /;
  FreeQ[{a,b},x] 

 ODEsecond[ a_. x_ Derivative[2][f_][x_]+(b_.+c_. x) f_[x_], f_[x_],x_,w_ ]:= 
   Module[ { newfun, answer,k=If[Znak[c],-Sqrt[-c/a],-Sqrt[c/a]] },
	answer = DSolveVictor[ E^(-x k)/x *
	(a x f''[x] + (b + c x) f[x]/.{
		f[x]->newfun[x] x E^(x k),
		f''[x]->D[newfun[x] x E^(x k),{x,2}]
		}) == 0, 
	newfun[x], x, DSolveConstants -> w ];
	If[ FreeQ[answer, FailODE],
	    {{ Rule[f[x], SSimplify[answer[[1,1,2]] E^(x k) x]] }},
	    FailODE
	]
   ] /;
 FreeQ[{a,c}, x]

 ODEsecond[ expr_, f_[x_], x_, w_ ] :=  ODEany[ expr,f[x],x,2,w ] 

 ODEany[ expr_, f_[x_], x_, order_ , w_ ] :=
   Block[ {Hypergeometric2F1,HypergeometricPFQ},
   Module[ {xnew,xold,r,answer},
	collectF = True;
	r = CheckCoef[ Expand[ CoefODE[
		CollectListDerivative[ Collect[expr, f[x]], f[x], order],
		f[x], order ]], x, order]; 
	If[ !FreeQ[ r, FailODE ], 
	    answer = FailODE,
	    If[ FreeQ[ r, ChangeVariable ],
		answer = ODEHyperType@@Join[ r, {f[x],x} ]; 
		If[ FreeQ[answer, FailODE],
		    answer = Sum[w@@{i} answer[[i]],{i,order}];
		    If[ FreeQ[answer,Integrate`sum] && order==2,
			answer = {{ f[x] -> ConstDelDSolve[answer,x,2,w]}},
			answer = {{ f[x] -> answer }}
		    ]
		],
		If[ !FreeQ[ r, ChangeVariable[] ],
		    r = ChangeVariables[ expr == 0,f,x,xnew,1/xnew ];
		    answer = DSolveVictor[PowerExpand[r],f[xnew],xnew,
				DSolveConstants -> w ]/.
				Rule[p_,q_]:>f[x]->q/. xnew->1/x,
		    If[ !FreeQ[ r, ChangeVariable[1] ],
			r = ChangeVariables[ expr == 0,f,x,xnew,xnew+1 ];
			answer = DSolveVictor[r,f[xnew],xnew,
				DSolveConstants -> w]/.
				Rule[p_,q_]:>f[x]->q/.xnew->x-1
			,
			{xold,xnew,val} = 
			  r/.ChangeVariable[ a_,b_,c_ ] :> {a,b,c};
			r = ChangeVariables[ expr == 0,f,x,xnew,val ]; 
			answer = 
			If[ !FreeQ[r, FailODE], FailODE,
			    DSolveVictor[PowerExpand[r],f[xnew],xnew,
				DSolveConstants->w]/.
				Rule[p_,q_]:>f[x]->q/.xnew->xold
			]
		   ]
		]
	    ]
	]; 
	answer
   ]] /; Depth[expr] < 10 && Length[expr] < 40

 ODEany[ __ ] := FailODE

 SSimplify[ expr_ ] := 
	Block[ {HypergeometricPFQ, HypergeometricU,
		Hypergeometric2F1, Hypergeometric1F1},
	Simplify[ expr ] 
	] /;
 And@@(FreeQ[Hold[expr],#]&/@{HypergeometricPFQ, HypergeometricU,
	Hypergeometric2F1, Hypergeometric1F1})

 SSimplify[ expr_ ] := Simplify[expr]	

 Attributes[ SSimplify ] = {HoldFirst}
	
 CollectListDerivative[ expr_, f_[x_], r_?Positive ] :=
	CollectListDerivative[ Collect[expr,f'[x]], f'[x], r-1 ]

 CollectListDerivative[ g_[expr__], _ , 0 ] := { expr } /; g===Plus 

 CollectListDerivative[ expr_, _ , 0 ] := { expr }

 CoefODE[ {a___, w_. Derivative[n_][f_][x_], b___}, f_[x_], n_ ] :=
	Prepend[ CoefODE[ {a,b}, f[x], n-1 ], w]

 CoefODE[ expr__,n_?Positive ] := Prepend[ CoefODE[ expr, n-1 ], 0]

 CoefODE[ {}, __ ] := {}

 CoefODE[ {a___, w_. f_[x_], b___}, f_[x_], 0 ] :=
	List[ If[ Length[ {a,b} ] == 0, w, FailODE ] ]

 CoefODE[ {}, __, -1 ] := { 0 }

 CoefODE[ {__}, __, -1 ] := { FailODE }
	
 FactorSolve[ a_ b_Plus, x_ ]:= If[ !FreeQ[b,x], b, a]

 FactorSolve[ __ ] :=  FailODE

 CheckCoef[ expr_, x_, n_ ] := {FailODE} /; !FreeQ[expr,FailODE]

 CheckCoef[ expr_, x_, n_ ] := {FailODE} /; 
 Length[ Flatten[ Cases[#,w_ x^k_/;Re[N[k]] < 0]&/@expr ] ] != 0

 CheckCoef[ expr_, x_, n_ ] := 
	( collectF=False; CheckCoef[Collect[expr,x], x, n] ) /;
 collectF

 CheckCoef[ {a_. x_^p_, rr__}, x_, n_ ] := 
	( olderdeg=False; ChangeVariable[] )/; 
 olderdeg && FreeQ[a,x] && Re[N[p-n]]>0

 CheckCoef[ {rr1___, a_. x_^p_. + b_. x_^q_. + c_, rr2___}, x_, n_ ] := 
   Module[ { root, var=Unique[var], oldvar ,const },
	Variable = False;
	Off[Solve::ifun, Solve::tdep];
	root = Solve[ a x^p + b x^q + c == var, x];
	On[Solve::ifun, Solve::tdep];
	If[ WordPresentQ[root], 
	    {FailODE},
	    ChangeVariable[ a x^p + b x^q + c, var, root[[1,1,2]] ]
	]
   ] /; Variable && FreeQ[{a,b,c},x] && Max[Re[{p,q}]]==2
   
 CheckCoef[ {rr___, a_. x_^p_. + b_. x_^q_.}, x_, n_ ] := 
   Module[ { root, varn, r },
	Variable = False;
	Off[Solve::ifun, Solve::tdep];
	root = Solve[ a x^p + b x^q == varn, x];
	On[Solve::ifun, Solve::tdep];
	If[ WordPresentQ[root], 
	    If[ (root=Factor[a x^p + b x^q]) =!= a x^p + b x^q,
		r = FactorSolve[root,x];
		If[ FreeQ[root,FailODE], 
		    Off[Solve::ifun, Solve::tdep];
		    root = Solve[ r == varn, x];
		    On[Solve::ifun, Solve::tdep];
		    If[ WordPresentQ[root], 
		        {FailODE},
		        ChangeVariable[ r, varn, root[[1,1,2]] ]
		    ],
		    {FailODE}
		],
		{FailODE}
	    ],
	    ChangeVariable[ a x^p + b x^q, varn, root[[1,1,2]] ]
	]
   ] /; Variable && FreeQ[{a,b,p,q},x]
	
 CheckCoef[ {rr1___, a_. x_^p_. + b_. x_^q_., rr2___}, x_, n_ ] := 
	If[ Max[p,q] < n - Length[{rr1}],
	    collectF = True;
	    CheckCoef[ Expand[x^(n-Max[p,q]-Length[{rr1}]) *
			{rr1,a x^p+b x^q,rr2}], x, n ],
	    If[ Min[p,q] =!= n - Length[{rr1}],
		collectF = True;
		CheckCoef[ Expand[x^(n-Min[p,q]-Length[{rr1}]) *
			  {rr1,a x^p+b x^q,rr2}], x, n 
		],
		If[ p == n - Length[{rr1}],
		    CheckCoefRest[ 
			Join[ CollectList[rr1,x], {{a,b}}, 
			      CollectList[rr2,x] ], 
		        x, n, Join[ {q-n+Length[{rr1}], n-Length[{rr1}], -b},
				    If[ Length[{rr1}] == 0, {a}, {} ] ]
		    ],
		    CheckCoefRest[ 
			Join[ CollectList[rr1,x], {{b,a}}, 
			      CollectList[rr2,x] ], 
		        x, n, Join[ {p-n+Length[{rr1}], n-Length[{rr1}], -a},
				    If[ Length[{rr1}] == 0, {b}, {} ] ]
		    ]
		]
	    ]
	] /;
 FreeQ[{a,b},x] && And@@((Im[N[#]]==0 && NumberQ[N[#]])&/@{p,q})

 CheckCoef[ {rr1___, a_ + b_. x_^q_., rr2___}, x_, n_ ] := 
   Module[ { answer },
	answer = 
	If[ n == Length[{rr1}],
	    CheckCoefRest[ 
		Join[ CollectList[rr1,x], {{a,b}}, CollectList[rr2,x] ],
		x, n, {q, 0, -b}
	    ], 
	    CheckCoef[Collect[Expand[x^(n-Length[{rr1}]) {rr1,a+b x^q,rr2}],x],
		x, n
	    ]
	];
	answer /; FreeQ[answer, FailODE]
   ] /;
 FreeQ[{a,b,q}, x]

 CheckCoef[ {rr___, a_. fun_ + b_.}, x_, n_ ] := 
   Module[ { varp,root },
        Variable = False;
	Off[Solve::ifun,Solve::tdep];
	root = Solve[ fun == varp, x ]; 
	On[Solve::ifun,Solve::tdep];
	root = 
	If[ WordPresentQ[root],
	    {FailODE},
	    ChangeVariable[ fun, varp, root[[1,1,2]] ]
	];
	root /; FreeQ[root, FailODE]
   ] /;Variable && n == Length[{rr}] && FreeQ[{a,b},x] && 
 !FreeQ[fun,x] && !MatchQ[fun,d_. x^_./;FreeQ[d,x]] && 
 !MatchQ[fun,d_. x^_. + t_/;FreeQ[{t,d},x]]

 CheckCoef[ {a_. x_^n_, rr___}, x_, n_ ] := 
	CheckCoefRest[ {{a,0}, rr}, x, n, {Fail,a} ] /;
 FreeQ[a,x]

 CheckCoef[ {a_. x_^k_., rr___}, x_, n_ ] := 
	CheckCoef[ x^(n-k) {a x^k, rr}, x, n ] /;
 FreeQ[a,x] && NumberQ[N[k]] && Im[N[k]]===0 

 CheckCoef[ {a_, rr___}, x_, n_ ] := 
	(Variable1 = False; CheckCoef[ Expand[x^n {a, rr}], x, n ] )/;
 Variable1 && FreeQ[a, x] && !FreeQ[{rr}, x]

 CheckCoef[ __ ] := {FailODE}

 CollectList[ list__, x_ ] := Collect[#,x]&/@{list}

 CollectList[ x_ ] := {}

 CheckCoefRest[ {rr1___, {a__}, rr2___}, x_, order_, l_ ] :=
	CheckCoefRest1[
	   Join[ CoefDevide[ {rr1}, x, order ],{{a}},
	         CoefDevide[ {rr2}, x, order - Length[{rr1}] - 1 ]
	   ],
	x, order, l
	]

 CoefDevide[ {}, x_, order_ ] := {}

 CoefDevide[ {w_, v___}, x_, 0 ] := {w}

 CoefDevide[ {w_, v___}, x_, order_ ] := 
	Prepend[ CoefDevide[{v},x,order-1], Collect[Expand[w x^(-order)],x] ]

 CheckCoefRest1[ {a_, rr___}, x_, n_, {c__} ] :=
	CheckCoefRest1[ {{a,0},rr}, x, n, {c,a} ] /;
 FreeQ[a,x] && Head[a] =!= List 

 CheckCoefRest1[ {a_. + b_. x_, rr___}, x_, n_, {1,c__} ] :=
	CheckCoefRest1[ {{a,b},rr}, x, n, {1,c,a} ] /;
 FreeQ[{a,b},x] 

 CheckCoefRest1[ {a_. + b_. x_^lambda_, rr___}, x_, n_, {lambda_,c__} ] :=
	CheckCoefRest1[ {{a,b},rr}, x, n, {lambda,c,a} ] /;
 FreeQ[{a,b},x]

 CheckCoefRest1[ {rr1__, a_. + b_. x_, rr2___}, x_, n_, {1,c__} ] :=
	CheckCoefRest1[ {rr1,{a,b},rr2}, x, n, {1,c} ] /;
 FreeQ[{a,b},x] 

 CheckCoefRest1[ {rr1__, a_. + b_. x_^l_, rr2___}, x_, n_, {l_,c__} ] :=
	CheckCoefRest1[ {rr1,{a,b},rr2}, x, n, {l,c} ] /;
 FreeQ[{a,b},x]

 CheckCoefRest1[ {rr1__, a_, rr2___}, x_, n_, c_ ] :=
	CheckCoefRest1[ {rr1,{a,0},rr2}, x, n, c ] /;
 FreeQ[a,x] && Head[a] =!= List

 CheckCoefRest1[ {rr1__, a_. + b_. x_^q_., rr2___}, x_, n_, {Fail, d_} ] :=
	CheckCoefRest1[ {rr1,{a,b},rr2}, x, n, {q,n-Length[{rr1}],-b,d} ] /;
 FreeQ[{a,b,q},x] 

 CheckCoefRest1[ a__, {Fail, w2_ } ] := {FailODE}

 CheckCoefRest1[ w__ ] := CheckCoefRest2[ w, {}, {} ]

 CheckCoefRest2[ { {a_, b_}, w___ }, x_, n_, l_, {up___}, {low___} ] :=
	CheckCoefRest2[ {w}, x, n, l, {up,a}, {low,b} ]

 CheckCoefRest2[ {w__}, x_, __ ] := (VariableC = False; ChangeVariable[1]) /;
 VariableC && Length[{w}]==2 && PolynomialQ[{w}[[1]],x]

 CheckCoefRest2[ {w__}, __ ] := {FailODE}

 CheckCoefRest2[ {}, x_, n_, {lambda_, p_, b_, 0}, up_, low_ ] := 
	ChangeVariable[]

 CheckCoefRest2[ {}, x_, n_, {lambda_, p_, b_, a_}, up_, low_ ] := 
	{ SolvePar[Reverse[up], n]/lambda, 
	  SolvePar[-Reverse[low], n]/lambda, x^lambda/a b lambda^(p-n), n }

 WordPresentQ[ word_ ] :=
 	Block[ { Solve },
		   !FreeQ[word,Roots] || !FreeQ[word,Solve] || 	
	            Length[word]==0
	]

 Attributes[WordPresentQ] = {HoldFirst}
 
 SolvePar[ {v_, list__}, n_ ] :=
   Module[ { x },
	Flatten[
	  Solve[
	    v + Sum[ {list}[[i]] Product[ x - j + 1, {j,i} ], {i,n} ] == 0, x
	  ] /. Rule[x, w_] :> w
	]
   ]

 ODEHyperType[ w__/; !FreeQ[{w},FailODE] ] := FailODE 

 ODEHyperType[ rootX_, __ ] := FailODE /; 
 Length[rootX] > 2 && LogCaseX[ rootX ]  

 ODEHyperType[ {x1_, x2_}, rootU_, psi_, 2, f_[x_], x_ ] :=
   Module[ { xx1=x1, xx2=x2, bb, dd },
	If[ Re[ N[x1-x2] ] < 0, xx1=x2; xx2=x1 ];
	If[ LogCase512[ xx1 - rootU ], 
	    { 
	      MeijerReduce[ {}, 1+rootU, {xx1}, {xx2}, 
				Exp[Pi I (Length[rootU]+1)] psi, x
	      ],
	      MeijerReduce[ {}, 1+rootU, {xx1,xx2}, {}, 
				Exp[Pi I Length[rootU]] psi, x
	      ] 
	    },
	    If[ LogCase512[ xx2 - rootU ], 
		{
		 If[ LogCase512[rootU-xx1],
		     MeijerReduce[ 1+rootU, {}, {xx1}, {xx2}, -psi, x ],
		     bb = Select[rootU,IntegerQ[xx1-#-1]&&Negative[xx1-#-1]&];
		     dd=Delete[rootU,Union[Flatten[Position[rootU,#]&/@bb,1]]];
		     MeijerReduce[ 1+dd, 1+bb, {xx1}, {xx2}, 
		    			Exp[Pi I (Length[bb]+1)] psi, x ]
		 ], 
		 MeijerReduce[ {}, 1+rootU, {xx1,xx2}, {}, 
				Exp[Pi I Length[rootU] ] psi, x
		 ] 
		},
		If[ LogCase512[1+rootU-xx2],
		    {
		    MeijerReduce[ 1+rootU, {}, {xx1}, {xx2}, -psi, x ], 
		    MeijerReduce[ 1+rootU, {}, {xx1,xx2}, {}, psi, x] 
		    }
		    ,
		    bb = Select[rootU, IntegerQ[xx2-#-1]&&Negative[xx2-#-1]& ];
		    dd =Delete[rootU,Union[Flatten[Position[rootU,#]&/@bb,1]]];
		    {
		    If[ LogCase512[1+rootU-xx1],
		        MeijerReduce[ 1+rootU, {}, {xx1}, {xx2},-psi, x ],
			MeijerReduce[ 1+dd, 1+bb, {xx1}, {xx2},  
		    			Exp[Pi I (Length[bb]+1)] psi, x]
		    ], 
		    MeijerReduce[ 1+dd, 1+bb, {xx1,xx2}, {}, 
		    			Exp[Pi I Length[bb]] psi, x] 
		    }
		]
	    ]
	]
   ] /;
 LogCaseX[ {x1,x2} ] 
 
 ODEHyperType[ rootX_, rootU_, psi_, order_, f_[x_], x_ ] :=
    Module[ { bb, dd },
	Table[ 
	     If[ LogCase512a[ rootX[[i]] - rootU ], 
		 bb = Select[rootU,IntegerQ[rootX[[i]]-#]&&
		     (Negative[rootX[[i]]-#] || SameQ[rootX[[i]]-#,0])& ];
		 dd =Delete[rootU,Union[Flatten[Position[rootU,#]&/@bb,1]]];
		 MeijerReduce[ 1+dd, 1+bb, {rootX[[i]]}, Drop[rootX,{i}], 
				Exp[Pi I (Length[bb]-1)] psi, x
		 ],
		 MeijerReduce[ 1+rootU, {}, {rootX[[i]]}, Drop[rootX,{i}], 
				Exp[Pi I] psi, x
		 ]
	     ],
	{i, order}
	] 
    ]

 MeijerReduce[ n_, p_, m_, q_, arg_, x_ ] :=
   Module[ { r1, r2 },
	r1 = ReducePar[ m,p ];
	r2 = ReducePar[ q,n ];
	MeijerPrelimCaseGlog[ r2[[2]], r1[[2]], r1[[1]], r2[[1]], arg ] /. 
		{ a_^k_ :> a^Expand[ k ] } //.
		{ (a_ x^b_.)^k_ :> a^k x^Expand[b k] }/.BesselComplexArg 
  ] 

 LogCaseX[ {x_, root__} ] := 
	If[ Or@@( IntegerQ[ # ]&/@( x-{root} )),
	    True,
	    LogCase[ {root} ]
	] 

 LogCaseX[ x_ ] := False

 LogCase512a[ {u1___, z_Integer, u2___} ] := True/;Negative[z]||z===0

 LogCase512a[ __ ] := False

 LogCase512[ {u1___, z_Integer?Positive, u2___} ] := False

 LogCase512[ __ ] := True

 ConstDelDSolve[ expr_?(!FreeQ[#,FailODE]&), __ ] := FailODE

 ConstDelDSolve[ expr_, x_ , 2, C_ ] := 
   Module[ { r = Expand[expr], r1, r2 }, 
	r1 = Cases[ Cases[ r,s_/; !FreeQ[s,C[1]] ]/C[1], s_/;FreeQ[s,C[_]] ];
	r2 = Cases[ Cases[ r,s_/; !FreeQ[s,C[2]] ]/C[2], s_/;FreeQ[s,C[_]] ];
	If[ Complement[r2,r1]==={},
	    r1 = Complement[r1,r2],
	    If[ Complement[r1,r2]==={},
		r2 = Complement[r2,r1] ]];
	C[1] SSimplify[Plus@@r1] + C[2] SSimplify[Plus@@r2] +
	Plus@@Cases[ r,s_/; FreeQ[s,C[1]] && FreeQ[s,C[2]]] 
   ] /; FreeQ[expr, Literal[Integrate[__]]]

ConstDelDSolve[ expr_, __ ] := expr

 MeijerPrelimCaseGlog[ n_, p_, m_, q_, arg_ ] :=
  Module[ { answer },
    answer = Integrate`MeijerPrelimCase[ n,p,m,q,arg ];
    If[ FreeQ[answer,Integrate`FailInt], 
        answer,
        answer = Integrate`MeijerPrelimCase[ 1-m,1-q,1-n,1-p,1/arg ];
        If[ FreeQ[answer,Integrate`FailInt], 
            answer,
            GfunToHyper[ n,p,m,q,arg ] 
        ]
    ]
  ]

 GfunToHyper[ n_, p_, {v1___, b_, v2___}, {v3___, a_, v4___}, arg_] :=
	(-1)^(a-b) GfunToHyper[ n,p,{a,v1,v2},{b,v3,v4},arg ] /;
 IntegerQ[Expand[a-b-1]] && NonNegative[Expand[a-b-1]]

 GfunToHyper[ n_, p_, m_, q_, arg_ ] := 
    If[ Length[n] + Length[p] == Length[m] + Length[q] && Length[n]!=0,
        GfunToHyper[ 1-m, 1-q, 1-n, 1-p, 1/arg],
	0
    ]/; Length[ m ] == 0
    
 GfunToHyper[ n_, p_, m_/;!LogarithmCase[ m ], q_, arg_ ] := 
	Sum[ FinalGfunToHyper[ m[[i]],n,p,Drop[ m,{i,i} ],q,arg],
            { i,Length[m] } 
	]   

 GfunToHyper[ n_, p_, m_, q_, arg_ ] := 
   Module[ { answer },
	answer = Integrate`MeijerLogCase[n,p,m,q,arg ];
	If[ FreeQ[answer,Integrate`FailInt],
	    answer/.{Integrate`sum :> Sum, Integrate`dummy :> Integrate`var},
	    FailODE
	]
   ]

 FinalGfunToHyper[ w_, n_, p_, m_, q_, arg_ ] :=
	arg^Together[w] *
	Module[ { ww },
	   SimplifyGamma[Limit[ SimplifyGamma[
			MultGamma[1 + ww - n] * MultGamma[m - ww] /
			(MultGamma[p - ww]  MultGamma[1 + ww - q]) ], 
		ww -> w]]
	]*
	HypergeometricPFQ[ Expand[1 + w - Join[n,p]],
			     Expand[1 + w - Join[m,q]],
			     arg (-1)^(Length[p]-Length[m]+1) ]

 LogarithmCase[ {___, v_ ,___, u_, ___} ] := True /; IntegerQ[Expand[u-v]]

 LogarithmCase[ {___} ] := False

 ReducePar[ {w1___, u_, w2___},{w3___, v_, w4___} ] := 
	ReducePar[ {w1, w2},{w3, w4} ] /; 
 u===v

 ReducePar[ w1_, w2_ ] := { w1, w2 }

 MultGamma[a_] := Times@@( Gamma[#]&/@Expand[a] ) 
 MultPochham[a_List,l_] := Times@@( Pochhammer[#,l]&/@Expand[a] ) 

 BesselComplexArg =
 {
   BesselJ[ v_, w_. Complex[0, a_] ] :> Exp[Pi v I/2] BesselI[v,a w],
   BesselI[ v_, w_. Complex[0, a_] ] :> Exp[-Pi v I/2] BesselJ[v,-a w]
 }

 ChangeVariables[ equation_, y_, x_, xnew_, expr_?(Less[Depth[#],9]&) ] :=
   equation/.
   {
   c_. Derivative[a_][y][x] :> 
   (c/.x->expr) DerivativeD[a, y, x, xnew, expr ],
   c_. y[x] :> (c/.x->expr) y[xnew]
   }/.
   Equal[ a_, b_] :> Equal[ Expand[a], Expand[b] ]/.InversTrig  

 ChangeVariables[ __ ] := FailODE
 
 DerivativeD[1, y_, x_, xnew_, expr_ ] := 
 Derivative[1][y][xnew]/D[expr,xnew]

 DerivativeD[ n_, y_, x_, xnew_, expr_ ] := 
 D[ DerivativeD[n-1,y,x,xnew,expr], xnew]/D[expr,xnew]

 InversTrig = {
  Cos[ArcSin[w_]]^m_. :> (1-w^2)^(m/2),
  Cos[ArcTan[w_]]^m_. :> 1/(1+w^2)^(m/2),
  Sec[ArcCos[w_]]^m_. :> w^(-m),
  Sec[ArcSin[w_]]^m_. :> (1-w^2)^(-m/2),
  Sec[ArcTan[w_]]^m_. :> (1+w^2)^(m/2),
  Sin[ArcCos[w_]]^m_. :> (1-w^2)^(m/2),
  Sin[ArcTan[w_]]^m_. :> w^m/(1+w^2)^(m/2),
  Csc[ArcSin[w_]]^m_. :> w^(-m),
  Csc[ArcCos[w_]]^m_. :> (1-w^2)^(-m/2),
  Csc[ArcTan[w_]]^m_. :> (1+w^2)^(m/2)/w^m,
  Tan[ArcSin[w_]]^m_. :> w^m/(1-w^2)^(m/2),
  Tan[ArcCos[w_]]^m_. :> (1-w^2)^(m/2)/w^m,
  Cot[ArcCos[w_]]^m_. :> w^m/(1-w^2)^(m/2),
  Cot[ArcSin[w_]]^m_. :> (1-w^2)^(m/2)/w^m,
  ArcTan[Tan[w_]]^m_. :> w^m
        }

 Kummer[ r1_ f_[x_] + r2_ Derivative[1][f_][x_] + 
	 r3_ Derivative[2][f_][x_], f_[x_], x_, w_ ] :=
   Module[ { a1, b1, a2, b2, a3, b3 },
     {a1,b1} = r1/.a_. + b_. x :> {a,b}/;FreeQ[{a,b},x];
     {a2,b2} = r2/.a_. + b_. x :> {a,b}/;FreeQ[{a,b},x];
     {a3,b3} = r3/.a_. + b_. x :> {a,b}/;FreeQ[{a,b},x];
   {{f[x]->E^(-(b2*x)/(2*b3) - (((b2^2 - 4*b1*b3)/b3^2)^(1/2)*(a3/b3 + x))/2 + 
      ((a3*b2 - a2*b3)*Log[a3 + b3*x])/(2*b3^2))*
   (a3/b3 + x)^((1 + ((a3*b2 - a2*b3 + b3^2)^2/b3^4)^(1/2))/2)*
   (w[1]*HypergeometricU[(-(a3*b2^2) + 2*a3*b1*b3 + a2*b2*b3 - 2*a1*b3^2)/
         (2*b3^3*((b2^2 - 4*b1*b3)/b3^2)^(1/2)) + 
        (1 + ((a3*b2 - a2*b3 + b3^2)^2/b3^4)^(1/2))/2, 
       1 + ((a3*b2 - a2*b3 + b3^2)^2/b3^4)^(1/2), 
       ((b2^2 - 4*b1*b3)/b3^2)^(1/2)*(a3/b3 + x)] + 
     w[2]*Hypergeometric1F1[(-(a3*b2^2) + 2*a3*b1*b3 + a2*b2*b3 - 
           2*a1*b3^2)/(2*b3^3*((b2^2 - 4*b1*b3)/b3^2)^(1/2)) + 
        (1 + ((a3*b2 - a2*b3 + b3^2)^2/b3^4)^(1/2))/2, 
       1 + ((a3*b2 - a2*b3 + b3^2)^2/b3^4)^(1/2), 
       ((b2^2 - 4*b1*b3)/b3^2)^(1/2)*(a3/b3 + x)]) }}/;b2^2 - 4*b1*b3=!=0
   ]/;
 MatchQ[r1, a_. + b_. x/;FreeQ[{a,b},x] ] &&
 MatchQ[r2, a_. + b_. x/;FreeQ[{a,b},x] ] &&
 MatchQ[r3, a_. + b_. x/;FreeQ[{a,b},x] ] 

 Kummer[ r1_ f_[x_] + a2_. Derivative[1][f_][x_] + 
	 r3_ Derivative[2][f_][x_], f_[x_], x_, w_ ] :=
   Module[ { a1, b1, a3, b3 },
     {a1,b1} = r1/.a_. + b_. x :> {a,b}/;FreeQ[{a,b},x];
     {a3,b3} = r3/.a_. + b_. x :> {a,b}/;FreeQ[{a,b},x];
  {{f[x] -> E^(-(((-4*b1)/b3)^(1/2)*(a3/b3 + x))/2 - 
         (a2*Log[a3 + b3*x])/(2*b3))*
      (a3/b3 + x)^((1 + ((-(a2*b3) + b3^2)^2/b3^4)^(1/2))/2)*
      (w[1]*HypergeometricU[(2*a3*b1*b3 - 2*a1*b3^2)/
            (2*((-4*b1)/b3)^(1/2)*b3^3) + 
           (1 + ((-(a2*b3) + b3^2)^2/b3^4)^(1/2))/2, 
          1 + ((-(a2*b3) + b3^2)^2/b3^4)^(1/2), 
          ((-4*b1)/b3)^(1/2)*(a3/b3 + x)] + 
        w[2]*Hypergeometric1F1[(2*a3*b1*b3 - 2*a1*b3^2)/
            (2*((-4*b1)/b3)^(1/2)*b3^3) + 
           (1 + ((-(a2*b3) + b3^2)^2/b3^4)^(1/2))/2, 
          1 + ((-(a2*b3) + b3^2)^2/b3^4)^(1/2), 
          ((-4*b1)/b3)^(1/2)*(a3/b3 + x)])}}
   ]/;
 MatchQ[r1, a_. + b_. x/;FreeQ[{a,b},x] ] &&
 MatchQ[r3, a_. + b_. x/;FreeQ[{a,b},x] ] && FreeQ[ a2, x]

 Kummer[ ___ ] := FailODE

 Unprotect[ HypergeometricU ]

 HypergeometricU[ a_, b_, x_ ] := 
   Simplify[Pi^(-1/2) Exp[x/2] x^(-a+1/2) BesselK[ a-1/2, x/2 ]] /;
 Expand[2 a - b] === 0

 Protect[ HypergeometricU ]

(*========================================================================*)

 SetAttributes[DSolve, { ReadProtected, Protected } ];

 End[ ] (* "Calculus`DSolve`Private` *)

 EndPackage[ ] (* "Calculus`DSolve`" *)

