# psolution is particular solution of linear differential equation;
# psolution2(b_,m_,c_,d_,y_,x_,n_):= if(m>0 and isfree(c/d,x), if(c==0, -int(b,x,m), -d/c));
# 	if(isfree(c1,x), psolution(a,n,b,m,c,c2,y,x,p)+psolution(c,c1,y,t,p), 
#	if(isfree(c2,x), psolution(a,n,b,m,c,c1,y,x,p)+psolution(c,c2,y,t,p), 



psolution2(a_,n_,b_,m_,c_,d_,y_,x_,p_):= if(has(d,t), psolution(c,d,y,t,p), if(has(d,x), psolution(b/a,m,c/a,d/a,y,x,n) ));
psolution(a_,n_,b_,m_,c_,c1_+c2_,y_,t_,p_):=
	if(has(c1,t) and hasnot(c1,x), if(has(c2,x) and hasnot(c2,t), psolution(b/a,m,c/a,c2/a,y,x,n),psolution(a,n,b,m,c,c2,y,t,p))+psolution(c,c1,y,t,p), 
	if(has(c2,t) and hasnot(c2,x), if(has(c1,x) and hasnot(c1,t), psolution(b/a,m,c/a,c1/a,y,x,n),psolution(a,n,b,m,c,c1,y,t,p))+psolution(c,c2,y,t,p), 
	if(has(c2,x) and hasnot(c2,t), psolution(a,n,b,m,c,c1,y,t,p)+psolution(b/a,m,c/a,c2/a,y,x,n), 
	psolution(a,n,b,m,c,c1,y,t,p)+psolution(a,n,b,m,c,c2,y,t,p) )));
psolution(a_,n_,b_,m_,c_,c1_*c2_,y_,t_,p_):= if(isfree(c1,x) and hasnot(c1,t), c1*psolution(a,n,b,m,c,c2,y,t,p), if(isfree(c2,x) and hasnot(c2,t), c2*psolution(a,n,b,m,c,c1,y,t,p), psolution(b,m,c,c1*c2,y,t,n) ));
psolution(a_,n_,b_,m_,c_,t*ax_,y_,t,p_):= if(p>1, psolution(b/a,m,c/a,ax/a,y,x,n)*t, if(c==0,psolution(a,n-1,b,t*ax,y,t,p) ));
psolution(a_,n_,b_,m_,c_,t*x_,y_,t,p_):=  if(p>1, psolution(b/a,m,c/a,x/a,y,x,n)*t,  if(c==0,psolution(a,n-1,b,t*x,y,t,p) ));
#psolution(a_,n_,b_,m_,c_,tt_*x,y_,t,p_):= if(p>1, psolution(b/a,m,c/a,x/a,y,x,n)*t, psolution(c,tt,y,t,p)*x+if(m==1,c*b*int(tt,t,p)/2,0) );
psolution(a_,n_,b_,m_,c_,tt_*x,y_,t,p_):= if(m>=1, psolution(-b,m,c,tt_*x,y,t,p) );
psolution(a_,n_,b_,m_,c_,exp(t)*xx_,y_,t,p_):= if(isconstant(a,b,c), psolution(b/a,m,(-1+c)/a,xx/a,y,x,n)*exp(t) );
psolution(a_,n_,b_,m_,c_,exp(t)*x,y_,t,p_):= if(isconstant(a,b,c), if(c==0 and m==1, (b+x)*exp(t), psolution(b,xx,y,x,1)*exp(t) )), 
psolution(a_,n_,b_,m_,c_,exp(t)*ax_,y_,t,p_):= if(isconstant(a,b,c), psolution(b/a,m,(-1+c)/a,ax/a,y,x,n)*exp(t) );
psolution(a_,n_,b_,m_,c_,exp(x)*tt_,y_,t,p_):= if(a-b-c==0 or isconstant(a,b,c), psolution(-a+b+c,tt,y,t,p)*exp(x) );
psolution(a_,n_,b_,m_,c_,exp(x)*at_,y_,t,p_):= if(a-b-c==0 or isconstant(a,b,c), psolution(-a+b+c,at,y,t,p)*exp(x) );

psolution(a_,n_,b_,m_,c_,c1_*c2_*c3_,y_,t_,p_):= if(isfree(c1,x) and hasnot(c1,t), c1*psolution(a,n,b,m,c,c2*c3,y,t,p), if(isfree(c2,x) and hasnot(c2,t), c2*psolution(a,n,b,m,c,c1*c3,y,t,p), psolution(b,m,c,c1*c2*c3,y,t,n) ));

psolution2(b_,m_,c_,d_,y_,x_,n_):=if(hasnot(d,y),if(isfree(c/d,x), if(c==0, -int(b,x,m),-d/c), if(b+c*x<>0 and hasnot(d/(b+c*x),x), -d*x/(b+c*x) )));
psolution(b_,1,c_,d_,y_,x_,2):=if(isfree(d/c,x), -d/c, if(b+c*x<>0 and hasnot(d/(b+c*x),x), -d*x/(b+c*x), if(isfree(b,c,y) and isfree(b,c,x), 
	if(b*b+4c==0, expand(exp(b/2*x)*x*int(exp(-b/2*x)*d,x)-exp(b/2*x)*int(d*exp(-b/2*x)*x,x)),
	expand((-exp((b-sqrt(b^2+4c))/2*x)*int(d*exp((-b+sqrt(b^2+4c))/2*x),x)+exp((b+sqrt(b^2+4c))/2*x)*int(d*exp((-b-sqrt(b^2+4c))/2*x),x))/sqrt(b*b+4c)) ))));

#psolution(b_,m_,c_,e^x_,y_,x_,n_):=if(isfree(b,c,y) and isfree(b,c,x), if(n==b*m, 1/n!*exp(x)*x^n, if(1-b-c==0, exp(x)*x/(n-b*m), exp(x)/(1-b-c) )));
#psolution(b_,m_,c_,exp(a_*x_),y_,x_,n_):=if(isfree(b,c,y) and isfree(b,c,x), if(a^n-b*a^m==c, exp(a*x)*x/(n*a^(n-1)-b*m*a^(m-1)), exp(a*x)/(a^n-b*a^m-c) ));
psolution(b_,m_,c_,exp(x_),y_,x_,n_):=if(isfree(b,c,y) and isfree(b,c,x), if(b+c==1, if(b*b== -4c and 2m==n, 1/m/n*exp(x)*x^2, exp(x)*x/(n-b*m) ), exp(x)/(1-b-c) ));
psolution(b_,m_,c_,exp(a_*x_),y_,x_,n_):=if(isfree(a,b,c,y) and isfree(a,b,c,x), if(a^n==b*a^m+c, if(b*b== -4c, a/n/m*exp(a*x)*x^2, exp(a*x)*x/(n*a^(n-1)-b*m*a^(m-1)) ), exp(a*x)/(a^n-b*a^m-c) ));
#psolution(b_,m_,c_,exp(x_+f_),y_,x_,n_):=if(isfree(b,c,f,y) and isfree(b,c,f,x), if(b+c==1, if(b*b== -4c, 1/n/m*exp(x+f)*x^2, exp(x+f)*x/(n-b*m) ), exp(x+f)/(1-b-c) ));
psolution(b_,m_,c_,exp(f_+x_),y_,x_,n_):=if(isfree(b,c,f,y) and isfree(b,c,x), if(b+c==1, if(b*b== -4c and 2m==n, 1/m/n*exp(x+f)*x^2, exp(x+f)*x/(n-b*m) ), exp(x+f)/(1-b-c) ));
psolution(b_,m_,c_,exp(f_+a_*x_),y_,x_,n_):=if(isfree(a,b,c,f,y) and isfree(a,b,c,f,x), if(a^n==b*a^m+c, if(b*b== -4c, a/n/m*exp(a*x+f)*x^2, exp(a*x+f)*x/(n*a^(n-1)-b*m*a^(m-1)) ), exp(a*x+f)/(a^n-b*a^m-c) ));
psolution(b_,m_,c_,sinh(x_),y_,x_,n_):=if(isfree(b,c,y) and isfree(b,c,x), psolution(b,m,c,exp(x)/2-exp(-x)/2,y,x,n));
psolution(b_,m_,c_,cosh(x_),y_,x_,n_):=if(isfree(b,c,y) and isfree(b,c,x), psolution(b,m,c,exp(x)/2+exp(-x)/2,y,x,n));
psolution(b_,m_,c_,sin(x_),y_,x_,n_):=if(isfree(b,c,y) and isfree(b,c,x), expand(totrig(psolution(b,m,c,i*exp(-i*x)/2-i*exp(i*x)/2,y,x,n))) );
psolution(b_,m_,c_,cos(x_),y_,x_,n_):=if(isfree(b,c,y) and isfree(b,c,x), expand(totrig(psolution(b,m,c,exp(i*x)/2+exp(-i*x)/2,y,x,n))) );
psolution(b_,m_,c_,sinh(k_*x_),y_,x_,n_):=if(isfree(b,c,y) and isfree(b,c,k,x), psolution(b,m,c,exp(k*x)/2-exp(-k*x)/2,y,x,n));
psolution(b_,m_,c_,cosh(k_*x_),y_,x_,n_):=if(isfree(b,c,y) and isfree(b,c,k,x), psolution(b,m,c,exp(k*x)/2+exp(-k*x)/2,y,x,n));
psolution(b_,m_,c_,sin(k_*x_),y_,x_,n_):=if(isfree(b,c,y) and isfree(b,c,k,x), expand(totrig(psolution(b,m,c,i*exp(-i*k*x)/2-i*exp(i*k*x)/2,y,x,n))) );
psolution(b_,m_,c_,cos(k_*x_),y_,x_,n_):=if(isfree(b,c,y) and isfree(b,c,k,x), expand(totrig(psolution(b,m,c,exp(i*k*x)/2+exp(-i*k*x)/2,y,x,n))) );
psolution(b_,m_,c_,sin(f_+x_),y_,x_,n_):=if(isfree(b,c,f,y) and isfree(b,c,f,x), expand(totrig(psolution(b,m,c,i*exp(f-i*x)/2-i*exp(i*x+f)/2,y,x,n))) );
psolution(b_,m_,c_,cos(f_+x_),y_,x_,n_):=if(isfree(b,c,f,y) and isfree(b,c,f,x), expand(totrig(psolution(b,m,c,exp(i*x+f)/2+exp(f-i*x)/2,y,x,n))) );
psolution(b_,m_,c_,sin(f_+k_*x_),y_,x_,n_):=if(isfree(b,c,k,f,y) and isfree(b,c,k,f,x), expand(totrig(psolution(b,m,c,i*exp(f-i*k*x)/2-i*exp(i*k*x+f)/2,y,x,n))) );
psolution(b_,m_,c_,cos(f_+k_*x_),y_,x_,n_):=if(isfree(b,c,k,f,y) and isfree(b,c,k,f,x), expand(totrig(psolution(b,m,c,exp(i*k*x+f)/2+exp(f-i*k*x)/2,y,x,n))) );
psolution(b_,n_,c_,x_,y_,x_,p_):=if(isfree(b,c,y) and isfree(b,c,x), if(n==1, -x/c+b/c^2, if(n>1, -x/c, b*x^n/n!/c^2-x/c-b^2/c^3 )));

psolution(a_,0,b_,d_,y_,t_,p_):= if(has(d,t), psolution(b,d,y,t,y), psolution(b/a,d/a,y,x,p) );
psolution(b_,1,0,d_,y_,x_,p_):= if(has(d,t), int(d,t,y), int(psolution(b,d,y,x,p-1),x) );



psolution(b_,m_,c_,c1_+c2_,y_,t_,n_):= 
	if(has(c1,t) and hasnot(c1,x), psolution(c,c1,y,t,n)+if(has(c2,x) and hasnot(c2,t), psolution(c/b,c2/b,n,x,m),psolution(b,m,c,c2,y,t,n) ), 
	if(has(c2,t) and hasnot(c2,x), psolution(c,c2,y,t,n)+if(has(c1,x) and hasnot(c1,t), psolution(c/b,c1/b,n,x,m),psolution(b,m,c,c1,y,t,n) ), 
	if(has(c2,x) and hasnot(c2,t), psolution(b,m,c,c1,y,t,n)+psolution(c/b,c2/b,n,x,m), 
	psolution(b,m,c,c1,y,t,n)+psolution(b,m,c,c2,y,t,n) )));
#psolution(b_,m_,c_,c1_+c2_,y_,x_,n_):= psolution(b,m,c,c1,y,x,n)+psolution(b,m,c,c2,y,x,n);
psolution(b_,m_,c_,c1_*c2_,y_,x_,n_):=if(isfree(c1,x), c1*psolution(b,m,c,c2,y,x,n), if(isfree(c2,x), c2*psolution(b,m,c,c1,y,x,n) ));

#psolution(b_,m_,c_,ax_*t,y_,t,p_):= if(isconstant(b,c), if(c==0, int(ax,x,m)*t/b-if(p==1,1/b/b*int(ax,x,2m),0),-1/2/c*x*exp(t)-1/4/c*exp(t)-exp(t)/2 ));
psolution(b_,m_,c_,ax_*t,y_,t,p_):= if(isconstant(b,c), if(c==0, int(ax,x,m)*t/b-if(p==1,1/b/b*int(ax,x,2m),0), if(p>1, psolution(c/b,ax/b,y,x,m)*t, psolution(c/b,ax/b,y,x,m)*t+c*int(ax/b,x,m)/2 )));
#psolution(b_,m_,c_,tt_*x,y_,t,p_):= if(isconstant(b,c), if(c==0, int(tt,t,p)*x-if(m==1,b*int(tt,t,2p),0), 1/b*exp(c*t) * integrate(exp(-c*t)* tt,t)*x+integrate(tt,t)/b/c/2 ));
#psolution(b_,m_,c_,tt_*x,y_,t,p_):= if(isconstant(b,c), if(c==0, int(tt,t,p)*x-if(m==1,b*int(tt,t,2p),0), if(m>1, psolution(c,tt,y,t,p)*x, -b*exp(c*t) * integrate(exp(-c*t)*tt,t)*x )));
psolution(b_,m_,c_,tt_*x,y_,t,p_):= if(isconstant(b,c), if(c==0, int(tt,t,p)*x-if(m>1,0,b*int(tt,t,2p)), if(m>1, psolution(c,tt,y,t,p)*x, psolution(c,tt,y,t,p)*x+c*int(tt,t,p)/2 )));
psolution(b_,m_,c_,t*x^n_,y_,t,p_):= if(isconstant(b,c), if(c==0, x^(m+n)*t/b/(n+m)-if(p>1,0,1/b/b*x^(n+2m)/(n+2m)/(n+m) )));
#psolution(b_,m_,c_,exp(t)*ax_,y_,t,p_):= if(isconstant(b,c), if(c==0, 1/b*int(exp(x/b)*ax,x,m)*exp(t+(-1/b)*x), int(ax,x,m)*exp(t) ));
#psolution(b_,m_,c_,exp(t)*xx_,y_,t,p_):= if(isconstant(b,c), if(c==0, 1/b*int(exp(x/b)*xx,x,m)*exp(t+(-1/b)*x), int(xx,x,m)*exp(t) ));
psolution(b_,m_,c_,exp(t)*ax_,y_,t,p_):= if(isconstant(b,c), psolution(-1/b+c/b,ax/b,y,x,m)*exp(t) );
psolution(b_,m_,c_,exp(t)*xx_,y_,t,p_):= if(isconstant(b,c), psolution(-1/b+c/b,xx/b,y,x,m)*exp(t) );
#psolution(b_,1,c_,exp(t)*x^n_,y_,t,p_):= if(isconstant(b,c), if(c==0, -(-b)^(1+n)*Gamma(1+n,(-1/b)*x)*exp(t+(-1/b)*x),-(-b/2)^(1+n)*Gamma(1+n,(-2/b)*x)*exp(t+(-2/b)*x) ));
#psolution(b_,m_,c_,exp(x)*at_,y_,t,p_):= if(isconstant(b,c), if(c==0, int(exp(b*t)*at,t,p)*exp((-b)*t+x), int(at,t,p)*exp(x) ));
#psolution(b_,m_,c_,exp(x)*tt_,y_,t,p_):= if(isconstant(b,c), if(c==0, int(exp(b*t)*tt,t,p)*exp((-b)*t+x), int(tt,t,p)*exp(x) ));
psolution(b_,m_,c_,exp(x)*at_,y_,t,p_):= if(isconstant(b,c) or c-b==0, psolution(-b+c,at,y,t,p)*exp(x) );
psolution(b_,m_,c_,exp(x)*tt_,y_,t,p_):= if(isconstant(b,c) or c-b==0, psolution(-b+c,tt,y,t,p)*exp(x) );
psolution(b_,m_,c_,exp(x)*t,y_,t,p_):= if(isconstant(b,c), if(c==0, 1/b*(if(p>1,0,if(p==1,-1/b,1/b-t^p/p!))+t)*exp(x),if(b==c, exp(x)*t*t/2,exp(x)*(-c*t+b*t-if(p>1,0,if(p==1,1,-p+t^p/p! )))/(b-c)^2 )));
psolution(b_,m_,c_,exp(t)*x,y_,t,p_):= if(isconstant(b,c), if(c==0, (if(m>1,0,if(m==1,-b,1/b-x^m/m!))+x)*exp(t),if(c==1,exp(t)/2*if(m>1,t*x,(2t*x-t*t)),exp(t)*((1-c)*x-if(m>1,0,if(m==1,b,-m+x^m/m! )))/(1-c)^2 )));
#psolution(b_,m_,c_,t*x,y_,t,p_):= if(isconstant(b,c), if(p>1, psolution(c/b,x/b,y,x,m)*t, if(m==1, if(c==0, int(t,t,p)*x-b*int(t,t,2p), -1/(c)*t*x-b/c/c*t-1/c/c*x+b*if(abs(c)==1,-2sgn(c),-sgn(c)/c/c) ), int(x,x,m)*t/b-1/b/b*int(x,x,2m) )));
psolution(b_,m_,c_,t*x,y_,t,p_):= if(isconstant(b,c), if(m>1, if(p>1, -1/c*t*x, psolution(c,t,y,t,p)*x), if(p>1, psolution(c,x,y,x,m)*t/b, 1/c*t*x+psolution(c,t,y,t,p)*x+psolution(c,x,y,x,m)*t/b-b*if(abs(c)==1,2sgn(c),sgn(c)/c/c)  )));
#psolution(b_,m_,0,t*x,y_,t,p_):= if(isconstant(b), if(p>1, if(m>1, 1/2*int(t,t,p)*x+1/2*int(x,x,m)*t/b,  psolution(0,x/b,y,x,m)*t), if(p<1,int(t,t,p)*x-b*int(t,t,2p),  if(m<1,int(x,x,m)*t/b-int(x,x,2m)/b^2 ))));
psolution(b_,m_,0,t*x,y_,t,p_):= if(isconstant(b), if(p>1, if(m>1, 1/2*int(t,t,p)*x+1/2*int(x,x,m)*t/b, int(x,x,m)*t/b-if(p>1,0,int(x,x,2m)/b^2 )),int(t,t,p)*x-if(m>1,0,b*int(t,t,2p)) ));
psolution(b_,m_,0,t*x,y_,t,1):= if(isconstant(b), if(m>1, b*t^2*x/2,int(x,x,m)*t/b-int(x,x,2m)/b^2 ));
psolution(b_,1,0,t*x,y_,t,p_):= if(isconstant(b), if(p>1, x*x*t/2/b,int(t,t,p)*x-b*int(t,t,2p) ));
psolution(b_,m_,c_,exp(t+x),y_,t,p_):= if(isconstant(b,c), exp(t+x)/(1+b-c) );

psolution(b_,m_,c_,c1_*c2_*c3_,y_,x_,n_):=if(isfree(c1,x), c1*psolution(b,m,c,c2*c3,y,x,n), if(isfree(c2,x), c2*psolution(b,m,c,c1*c3,y,x,n) ));

psolution(b_,1,c_,sin(x_),y_,x_,2):=if(isfree(b,c,y) and isfree(b,c,x), - (c+1)*sin(x)/(b^2 + (c + 1)^2) + b*cos(x)/(b^2 + (c + 1)^2) );
psolution(b_,1,c_,cos(x_),y_,x_,2):=if(isfree(b,c,y) and isfree(b,c,x), - b*sin(x)/(b^2 + (c + 1)^2) - (c+1)*cos(x)/(b^2 + (c + 1)^2) );
psolution(b_,1,c_,sin(k_*x_),y_,x_,2):=if(b==0 and c== -k*k, sgn(k)*cos(k*x)*x/c, if(isfree(b,c,k,y) and isfree(b,c,k,x), b*k*cos(k*x)/(b^2*k^2 + c^2 + 2c*k^2+k^4) - (c+k^2)*sin(k*x)/(b^2*k^2 + c^2 + 2c*k^2+k^4)  ));
psolution(b_,1,c_,cos(k_*x_),y_,x_,2):=if(b==0 and c== -k*k, -sin(k*x)*x/c, if(isfree(b,c,k,y) and isfree(b,c,k,x), -b*k*sin(k*x)/(b^2*k^2 + c^2 + 2c*k^2+k^4) - (c+k^2)*cos(k*x)/(b^2*k^2 + c^2 + 2c*k^2+k^4)  ));
psolution(b_,1,c_,sin(k_+x_),y_,x_,2):=if(isfree(b,c,k,y) and isfree(b,c,k,x), - (c+1)*sin(x+k)/(b^2 + (c + 1)^2) + b*cos(x+k)/(b^2 + (c + 1)^2) );
psolution(b_,1,c_,cos(k_+x_),y_,x_,2):=if(isfree(b,c,k,y) and isfree(b,c,k,x), - b*sin(x+k)/(b^2 + (c + 1)^2) - (c+1)*cos(x+k)/(b^2 + (c + 1)^2) );
psolution(b_,1,c_,sin(f_+k_*x_),y_,x_,2):=if(b==0 and c== -k*k, sgn(k)*cos(f+k*x)*x/c, if(isfree(b,c,k,f,y) and isfree(b,c,k,f,x),  b*k*cos(k*x+f)/(b^2*k^2 + c^2 + 2c*k^2+k^4) - (c+k^2)*sin(k*x+f)/(b^2*k^2 + c^2 + 2c*k^2+k^4)  ));
psolution(b_,1,c_,cos(f_+k_*x_),y_,x_,2):=if(b==0 and c== -k*k, -sin(f+k*x)*x/c, if(isfree(b,c,k,f,y) and isfree(b,c,k,f,x), -b*k*sin(k*x+f)/(b^2*k^2 + c^2 + 2c*k^2+k^4) - (c+k^2)*cos(k*x+f)/(b^2*k^2 + c^2 + 2c*k^2+k^4)  ));
psolution(b_/x_,1,x_^(-2),d_,y_,x_,2):=if(isfree(d,y), int(replace(d*x^2,x,exp(x)),x) );
psolution(1/x_,1,c_*x_^(-2),d_,y_,x_,2):=if(isfree(d,y), int(replace(d*x^2,x,exp(x)),x) );
psolution(b_/x_,1,c_*x_^(-2),d_,y_,x_,2):=if(isfree(d,y), int(replace(d*x^2,x,exp(x)),x) );
#psolution(b_,1,c_*x_^(-2),d_*x_^(-2),y_,x_,2):=if(isfree(c,d,x), -d/c);
#psolution(b_,1,c_*x_^(-2),x_^(-2),y_,x_,2):=if(isfree(c,x), -1/c);
#psolution(b_,1,x_^(-2),d_*x_^(-2),y_,x_,2):=if(isfree(d,x), -d);
psolution(b_,1,x_,f_,y_,x_,2):= if(isconstant(b), Bi(x+1/4b*b)*exp(b*x/2)*int(Ai(x+1/4b*b)*exp(-b*x/2)*f,x)-Ai(x+1/4b*b)*exp(b*x/2)*int(Bi(x+1/4b*b)*exp(-b*x/2)*f,x));
psolution(b_,1,c_*x_,f_,y_,x_,2):= if(isconstant(b,c), Bi((c*x+1/4b*b)/b^(2/3))*exp(b*x/2)*int(Ai((c*x+1/4b*b)/b^(2/3))*exp(-b*x/2)*f,x)/cbrt(c)-Ai((c*x+1/4b*b)/b^(2/3))*exp(b*x/2)*int(Bi((c*x+1/4b*b)/b^(2/3))*exp(-b*x/2)*f,x)/cbrt(c));

psolution(1,1,x_,f_,y_,x_,2):= Bi(x+1/4)*exp(x/2)*int(Ai(x+1/4)*exp(-x/2)*f,x)-Ai(x+1/4)*exp(x/2)*int(Bi(x+1/4)*exp(-x/2)*f,x);
psolution(1,1,c_*x_,f_,y_,x_,2):= if(isconstant(c), Bi((b*x+1/4)/b^(2/3))*exp(x/2)*int(Ai((b*x+1/4)/b^(2/3))*exp(-x/2)*f,x)/cbrt(c)-Ai((b*x+1/4)/b^(2/3))*exp(x/2)*int(Bi((b*x+1/4)/b^(2/3))*exp(-x/2)*f,x)/cbrt(c));
psolution(x_,1,-x_,exp(x_),y_,x_,2):=exp(x);
psolution(-x_,1,x_,exp(x_),y_,x_,2):=exp(x);
psolution(1,1,0,sin(x_),y,x_,2):= cos(x)/2-sin(x)/2;
psolution(1,1,0,cos(x_),y,x_,2):= -cos(x)/2-sin(x)/2;
psolution(1,1,0,exp(x_),y,x_,2):= exp(x)*x;
psolution(1,1,0,exp(b_*x_),y,x_,2):= if(isconstant(b), exp(b*x)/b );
psolution(b_,1,0,x_,y,x_,2):= if(isconstant(b), -1/b/2*x*x-x/b^2 );
psolution(b_,1,0,x_^2,y,x_,2):= if(isconstant(b), -2x/b^3-x^2/b^2-x^3/b/3 );

psolution(b_,m_,c_,c_,y_,x_,n_):= if(m>0, -1);
psolution(b_,m_,c_,bb_*c_,y_,x_,n_):=if(isfree(bb,c,y) and isfree(bb,x),  -bb);
psolution(b_,m_,a_*c_,c_,y_,x_,n_):=if(isfree(a,c,y) and isfree(a,x),  -1/a);
psolution(b_,m_,a_*c_,bb_*c_,y_,x_,n_):=if(isfree(a,bb,c,y) and isfree(a,bb,x), -bb/a);


#psolution2(c_,d_,y_,x_,n_):=if(isfree(c,d,y), if(isfree(c,d,x), -d/c, if(n==2, if(isfree(c,x), 
	if(c<0, expand((cos(sqrt(-c)*x)*int(d*sin(sqrt(-c)*x),x)-sin(sqrt(-c)*x)*int(d*cos(sqrt(-c)*x),x))/2/sqrt(-c)),
		expand((exp(sqrt(c)*x)*int(d*exp(-sqrt(c)*x),x)-exp(-sqrt(c)*x)*int(d*exp(sqrt(c)*x),x))/2/sqrt(c)) )),
	if(isinteger(n), block(dsolve:=integrate(c,x), expand(exp(dsolve)*integrate(d*exp(-dsolve),x))),
	 int((x-z)^(n-1)*replace(d,x,z)*mittag(n,n,replace(int(c,x,n),x,x-z)*n!), z,0,x) ))));
psolution2(c_,d_,y_,x_,n_):=if(isfree(c,d,y), if(isfree(c,d,x), -d/c, if(n==2, if(isfree(c,x), 
	if(c<0, expand((cos(sqrt(-c)*x)*int(d*sin(sqrt(-c)*x),x)-sin(sqrt(-c)*x)*int(d*cos(sqrt(-c)*x),x))/2/sqrt(-c)),
		expand((exp(sqrt(c)*x)*int(d*exp(-sqrt(c)*x),x)-exp(-sqrt(c)*x)*int(d*exp(sqrt(c)*x),x))/2/sqrt(c)) )),
	block(dsolve:=integrate(c,x,n)*n!, expand(integrate(d/E(n,dsolve),x,n)*E(n,dsolve))) )));
#psolution2(c_,d_,y_,x_,n_):=if(isfree(c,d,y) and n<1, block(gs:=gsolution(c,y,x,n)/C_1, int(simplify(d/gs),x,n)*gsolution(c,y,x,n)/C_1) );
#psolution2(c_,d_,y_,x_,n_):=if(isfree(c,d,y) and n<=1, block(gs:=gsolution(c,y,x,n)/C_1, int(simplify(d/gs),x,n)*gs) );
#psolution(c_,d_,y_,x_,2):=if(isfree(c,d,y) and isfree(c,x),if(c<0, cos(sqrt(-c)*x)*int(d*sin(sqrt(-c)*x),x)/2/sqrt(-c)-sin(sqrt(-c)*x)*int(d*cos(sqrt(-c)*x),x)/2/sqrt(-c),
	expand(exp(sqrt(c)*x)*int(d*exp(-sqrt(c)*x),x)/2/sqrt(c)-exp(-sqrt(c)*x)*int(d*exp(sqrt(c)*x),x)/2/sqrt(c)) ));

psolution(c_,g_,y_,x_,1):=if(hasnot(g/c,x), -g/c, if(isfree(c,g,y), block(dsolve:=integrate(c,x), exp(dsolve)*(integrate(g*exp(-dsolve),x))) ));

psolution(c_,c1_+c2_,yy_,x_,n_):=psolution(c,c1,yy,x,n)+psolution(c,c2,yy,x,n);
#psolution(c_,c1_+c2_,yy_,x_,n_):=if(isconstant(c,c1), psolution(c,c2,yy,x,n)-c1/c, if(isfree(c1,x), psolution(c,c2,yy,x,n)+if(c==0, int(c1,t,yy),psolution(c,c1,y,t,yy)), if(isfree(c2,x), psolution(c,c1,yy,x,n)+if(c==0, int(c2,t,yy),psolution(c,c2,y,t,yy)), psolution(c,c1,yy,x,n)+psolution(c,c2,yy,x,n) )));
psolution(c_,c1_*c2_,y_,x_,n_):=if(isfree(c1,x), c1*psolution(c,c2,y,x,n), if(isfree(c2,x), c2*psolution(c,c1,y,x,n) ));
psolution(b_,c_*d_,1,x_,1):= if(hasnot(c,x) and has(d,x) and hasnot(d,t), expand(exp(b*x)*int(exp(-b*z)*replace(c,t,t+x-z)*replace(d,x,z),z,0,x)), if(has(d,t) and hasnot(d,x) and has(c,x) and hasnot(c,t), expand(exp(b*x)*int(exp(-b*z)*replace(d,t,t+x-z)*replace(c,x,z),z,0,x)) ));
psolution(c_,c1_*c2_*c3_,y_,x_,n_):=if(isfree(c1,x), c1*psolution(c,c2*c3,y,x,n), if(isfree(c2,x), c2*psolution(c,c1*c3,y,x,n) ));

psolution(c_,E(a_,x_^a_),y_,x_,n_):=if(isfree(c,y) and isfree(c,x) and a==n, if(c==1, 2/sqrt(pi)*E(a,x^a_)*x^a, E(a,x^a_)/(1-c) ));
psolution(c_,E(a_,b_*x_^a_),y_,x_,n_):=if(isfree(b,c,y) and isfree(b,c,x) and a==n, if(c==b, 2/sqrt(pi)*E(a,b*x^a)*x^a,E(a,b*x^a)/(b-c) ));

psolution(c_,sinh(x_),y_,x_,n_):=if(isfree(c,y) and isfree(c,x), psolution(c,exp(x)/2-exp(-x)/2,y,x,n));
psolution(c_,cosh(x_),y_,x_,n_):=if(isfree(c,y) and isfree(c,x), psolution(c,exp(x)/2+exp(-x)/2,y,x,n));
psolution(c_,sin(x_),y_,x_,n_):=if(isfree(c,y) and isfree(c,x) and n<>1, expand(totrig(psolution(c,i*exp(-i*x)/2-i*exp(i*x)/2,y,x,n))) );
psolution(c_,cos(x_),y_,x_,n_):=if(isfree(c,y) and isfree(c,x) and n<>1, expand(totrig(psolution(c,exp(i*x)/2+exp(-i*x)/2,y,x,n))) );
psolution(c_,sinh(k_*x_),y_,x_,n_):=if(isfree(c,y) and isfree(c,x), psolution(c,exp(k*x)/2-exp(-k*x)/2,y,x,n));
psolution(c_,cosh(k_*x_),y_,x_,n_):=if(isfree(c,y) and isfree(c,x), psolution(c,exp(k*x)/2+exp(-k*x)/2,y,x,n));
psolution(c_,sin(k_*x_),y_,x_,n_):=if(isfree(c,y) and isfree(c,x) and n<>1, expand(totrig(psolution(c,i*exp(-i*k*x)/2-i*exp(i*k*x)/2,y,x,n))) );
psolution(c_,cos(k_*x_),y_,x_,n_):=if(isfree(c,y) and isfree(c,x) and n<>1, expand(totrig(psolution(c,exp(i*k*x)/2+exp(-i*k*x)/2,y,x,n))) );
psolution(c_,sin(f_+x_),y_,x_,n_):=if(isfree(c,y) and isfree(c,x) and n<>1, expand(totrig(psolution(c,i*exp(-i*x-f*i)/2-i*exp(i*x+f*i)/2,y,x,n))) );
psolution(c_,cos(f_+x_),y_,x_,n_):=if(isfree(c,y) and isfree(c,x) and n<>1, expand(totrig(psolution(c,exp(i*x+f*i)/2+exp(-i*x-f*i)/2,y,x,n))) );
psolution(c_,sin(f_+k_*x_),y_,x_,n_):=if(isfree(c,y) and isfree(c,x) and n<>1, expand(totrig(psolution(c,i*exp(-i*k*x-f*i)/2-i*exp(i*k*x+f*i)/2,y,x,n))) );
psolution(c_,cos(f_+k_*x_),y_,x_,n_):=if(isfree(c,y) and isfree(c,x) and n<>1, expand(totrig(psolution(c,exp(i*k*x+f*i)/2+exp(-i*k*x-f*i)/2,y,x,n))) );
psolution(c_,exp(a_*x_),y_,x_,n_):=if(isfree(a,c,y) and isfree(c,x), if(a^n==c, exp(a*x)*x/n/a^(n-1), exp(a*x)/(a^n-c) ));
psolution(c_,exp(a_+x_),y_,x_,n_):=if(isfree(a,c,y) and isfree(a,c,x), if(c==1, exp(a+x)*x/n, exp(a+x)/(1-c) ));
psolution(c_,exp(b_+a_*x_),y_,x_,n_):=if(isfree(b,c,y) and isfree(b,c,x), if(a^n==c, exp(a*x+b)*x/n/a^(n-1), exp(a*x+b)/(a^n-c) ));
psolution(c_,exp(x_)*x_,y_,x_,n_):=if(isfree(c,y) and isfree(c,x), if(c==1, exp(x)*x^2/2/n-(n-1)*exp(x)*x/2/n, -n*exp(x)/(c-1)^2-exp(x)*x/(c-1) ));
psolution(c_,exp(a_*x_)*x_,y_,x_,n_):=if(isfree(a,c,y) and isfree(a,c,x), if(a^n==c, exp(a*x)*x^2/(2n*a^(n-1))-(n-1)*exp(a*x)*x/(2n*a^n), -a^(n-1)*n*exp(a*x)/(c-a^n)^2-exp(a*x)*x/(c-a^n) ));
psolution(c_,exp(d_+x_)*x_,y_,x_,n_):=if(isfree(c,d,y) and isfree(c,d,x), if(c==1, exp(d+x)*x^2/2/n-(n-1)*exp(d+x)*x/2/n, -n*exp(d+x)/(c-1)^2-exp(d+x)*x/(c-1) ));
psolution(c_,exp(b_+a_*x_)*x_,y_,x_,n_):=if(isfree(a,b,c,y) and isfree(a,b,c,x), if(a^n==c, exp(a*x+b)*x^2/(2n*a^(n-1))-(n-1)*exp(a*x+b)*x/(2n*a^n), -a^(n-1)*n*exp(a*x+b)/(c-a^n)^2-exp(a*x+b)*x/(c-a^n) ));
psolution(c_,exp(x_),y_,x_,n_):=if(isfree(c,y) and isfree(c,x), if(c==1, exp(x)*x/n, exp(x)/(1-c) ));

#psolution(b_,x_^n_, y_,x_,q_) := If(isfree(b,x) and isfree(b,y), If(n==q, -x^q/b-q!/b^2, If(q>0, n!*x^(n+q)*mittag(q,1+n+q,b*x^q), Gamma(n+1)*x^n*mittag(-q,n+1,1/b*x^(-q)) )));
psolution(b_,x_^m_, y_,x_,q_) := If(isfree(b,x) and isfree(b,y) and m>0, if(q==m, -x^m/b-m!, if(n(remainder(m,q))==0, -sum(b^(-1-k/q)*m!/(m-k)!*x^(m-k),k,0,m,q) )));
psolution(c_,x_,y_,x_,n_):=if(isfree(c,y) and isfree(c,x), if(n==1, -x/c-1/c^2, if(n== -1,1,if(n<0, 1+x^(1+n)/(1+n)!,if(n>1, -x/c, -x/c-x^n/n!/c^2-1/c^3 )))));

psolution(x_,f_,y_,x_,2):= Bi(x)*int(Ai(x)*f,x)-Ai(x)*int(Bi(x)*f,x);
psolution(c_*x_,f_,y_,x_,2):= if(isconstant(c),Bi(cbrt(c)*x)*int(Ai(cbrt(c)*x)*f,x)/cbrt(c)-Ai(cbrt(c)*x)*int(Bi(cbrt(c)*x)*f,x)/cbrt(c));
#psolution(-1,d_,y_,x_,2):= expand(int(d*cos(x),x)*sin(x)-int(d*sin(x),x)*cos(x));
psolution(c_,sin(x_),y_,x_,2):=if(isfree(c,x), if(c== -1, -cos(x)*x/2, -sin(x)/(c + 1) ));
psolution(c_,cos(x_),y_,x_,2):=if(isfree(c,x), if(c== -1,  sin(x)*x/2, -cos(x)/(c + 1) ));
psolution(c_,tan(x_),y_,x_,2):=if(isfree(c,x), if(c== -1, -cos(x)*log(sec(x)+tan(x)) ));
psolution(c_,cot(x_),y_,x_,2):=if(isfree(c,x), if(c== -1,  sin(x)*log(sec(x)+tan(x)) ));
psolution(c_,sin(d_*x_),y_,x_,2):=if(c==-d*d, sgn(d)*cos(d*x)*x/c, if(isfree(c,d,x), -sin(d*x)/(c+d^2) ));
psolution(c_,cos(d_*x_),y_,x_,2):=if(c==-d*d, -sin(d*x)*x/c, if(isfree(c,d,x), -cos(d*x)/(c+d^2) ));
psolution(c_,sin(b_+x_),y_,x_,2):=if(isfree(b,c,x), if(c== -1, -cos(x+b)*x/2, -sin(x+b)/(c + 1) ));
psolution(c_,cos(b_+x_),y_,x_,2):=if(isfree(b,c,x), if(c== -1,  sin(x+b)*x/2, -cos(x+b)/(c + 1) ));
psolution(c_,sin(f_+d_*x_),y_,x_,2):=if(c==-d*d, sgn(d)*cos(f+d*x)*x/c, if(isfree(c,d,f,x), -sin(d*x+f)/(c+d^2) ));
psolution(c_,cos(f_+d_*x_),y_,x_,2):=if(c==-d*d, -sin(f+d*x)*x/c, if(isfree(c,d,f,x), -cos(d*x+f)/(c+d^2) ));


psolution(c_,c_,y_,x_,n_):=if(isfree(c,y), -1);
psolution(c_,b_*c_,y_,x_,n_):=if(isfree(b,c,y) and isfree(b,x),  -b);
psolution(a_*c_,c_,y_,x_,n_):=if(isfree(a,c,y) and isfree(a,x),  -1/a);
psolution(a_*c_,b_*c_,y_,x_,n_):=if(isfree(a,b,c,y) and isfree(a,b,x), -b/a);
psolution(0,c_,y_,x_,n_):=if(isfree(c,y), int(c,x,n) );

psolution(b_/x_,c_,y_,x_,1):= int(c/x^b,x)*x^b;
psolution(1/x_,c_,y_,x_,1):= int(c/x,x)*x;