read "bessellist";
read "equiv";

transfo:=proc(L,a)  
        local i,f;  
        global x,Dx;
        f:=add(mult( subs(x=a,coeff(L,Dx,i)) , (1/diff(a,x) * Dx)$i
                ),i=0..degree(L,Dx));
        sort(collect(f/lcoeff(f,Dx),Dx,factor),Dx)
end:

SIMP := proc(a)
	collect(a, indets(a,function), factor)
end:



solvebessel:=proc(L)
	
	#option trace;
	local Lb,list,i,niu,L0,yb,y0,res,G,y;
	global x, Dx;
	res:={};
	Lb:= x^2*Dx^2+x*Dx-(x^2+niu^2);
	yb:=diffop2de(Lb,y(x));
	yb:=rhs(dsolve(yb));
	
	list:=bessellist(L);
	if list <> {} then
		for i in list do
			L0:=subs(niu=i[2],Lb);
			y0:= subs(niu=i[2],yb);
			L0:= transfo(L0,i[1]);
			y0:=subs(x=i[1],y0);
			G:=equiv(L0,L);
			if G<>0 then 
				res := eval( diffop2de(G,y(x)), y(x) = y0);
				res := collect(res,{_C1,_C2}, SIMP);
				break;
			fi;
		od;
	fi;
	res;
end: