# 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));


psolution2(a_,n_,b_,m_,c_,d_,y_,x_,p_):= if(has(d,t), psolution(b,m,c,d,y,t,p), if(has(d,x), psolution(b/a,m,c/a,d/a,y,x,n) ));
psolution(a_,n_,b_,1,0,d_,y_,tx,p_):= if(has(d,t), int(d,t,y), int(psolution(b,d,y,x,n-1),x) );
psolution(a_,n_,b_,m_,c_,c1_*c2_,y_,tx_,p_):= if(isfree(c1,x) and hasnot(c1,t), c1*psolution(a,n,b,m,c,c2,y,tx,p), if(isfree(c2,x) and hasnot(c2,t), c2*psolution(a,n,b,m,c,c1,y,tx,p) ));

psolution(a_,n_,b_,m_,c_,t*ax_,y_,tx,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_,tx,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_,tx,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_,tx,p_):= if(m>=1, psolution(-b,m,c,tt*x,y,t,p) );
psolution(a_,n_,b_,m_,c_,exp(x)*tt_,y_,tx,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_,tx,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_,exp(t)*xx_,y_,tx,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)*ax_,y_,tx,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(t)*x,y_,tx,p_):= if(isconstant(a,b,c), if(c==0 and m==1, (b+x)*exp(t), psolution(c-1,x,y,x,m)*exp(t) ));
#psolution(a_,n_,b_,m_,c_,exp(t)*x,y_,tx,p_):= if(isconstant(a,b,c), if(c==0 and m==1, (b+x)*exp(t), if(c==1, psolution(b,m,0,x,y,x,n)*exp(t), if(m>=1, psolution(-b,m,c,exp(t)*x,y,t,p) ))));
psolution(a_,n_,b_,m_,c_,exp(t+x),y_,tx,p_):= if(isconstant(b,c), exp(t+x)/(1+a-b-c) );

psolution(a_,n_,b_,m_,c_,c1_*c2_*c3_,y_,t_,p_):= if(isconstant(c1), c1*psolution(a,n,b,m,c,c2*c3,y,t,p), if(isconstant(c2), c2*psolution(a,n,b,m,c,c1*c3,y,t,p), if(isconstant(c3), c3*psolution(b,m,c,c1*c2,y,t,n) )));
psolution(a_,n_,b_,m_,c_,c1_+c2_,y_,x_,p_):=psolution(a,n,b,m,c,c1,y,x,p)+psolution(a,n,b,m,c,c2,y,x,p);
psolution(a_,n_,b_,m_,c_,c1_+c2_,y_,tx,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,tx,p))+psolution(c,c1,y,t,p), 
	if(has(c1,x) and hasnot(c1,t), psolution(a,n,b,m,c,c2,y,tx,p)+psolution(b/a,m,c/a,c1/a,y,x,n),
	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,tx,p))+psolution(c,c2,y,t,p), 
	if(has(c2,x) and hasnot(c2,t), psolution(a,n,b,m,c,c1,y,tx,p)+psolution(b/a,m,c/a,c2/a,y,x,n), 
	psolution(a,n,b,m,c,c1,y,tx,p)+psolution(a,n,b,m,c,c2,y,tx,p)
	))));



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) )));
psolution2(b_,1,c_,d_,y_,x_,2):=if(isfree(c/d,x), if(c==0,-b*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_,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_,m_,c_,exp(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_,x_^n_,y_,x_,p_):=if(c<>0 and isfree(b,c,y) and isfree(c,x), if(m==n, -x^n/c+b/c^2, if(m>n, -x^n/c, b*x^n/n!/c^2-x/c-b^2/c^3 )));
psolution(b_,m_,c_,x_,y_,x_,p_):=if(c<>0 and isfree(b,c,y) and isfree(c,x), if(m==1, -x/c+b/c^2, if(m>1, -x/c, b*x^m/m!/c^2-x/c-b^2/c^3 )));
psolution(b_,mm_,x_^m_,x_^n_,y_,x_,p_):=if(mm>n-m and n-m>0, -x^(n-m));
psolution(b_,mm_,x_,x_^n_,y_,x_,p_):=if(mm>n-1 and n>1, -x^(n-1));

psolution(b_,m_,0,d_,y_,tx,m_):= if(isconstant(b), if(has(d,t-x) or has(d,x-t), if(m==1, if(b==1, d*x, int(d,x,m)/(b-1)), int(d,x,m)/(1+b)), if(has(d,t+x), if(b== -1, d*x, int(d,x,m)/(1+b)  ))));
 
#psolution(b_,m_,c_,c1_*c2_,y_,t_,n_):=if(hasnot(c1,t) and hasnot(c1,x), c1*psolution(b,m,c,c2,y,t,n), if(hasnot(c2,t) and hasnot(c2,x), c2*psolution(b,m,c,c1,y,t,n),
	if(m==1 and n==1 and has(c1,t) and hasnot(c1,x), if(c==0, int(replace(c1,t,t-x/b+z/b)*replace(c2,x,z),z,1,x)/b,
	int(exp(-c*z/b)*replace(c1,t,t-x/b+z/b)*replace(c2,x,z),z,1,x)*exp(c*x/b)/b ))));
psolution(b_,m_,c_,c1_*c2_,y_,t_,n_):=if(hasnot(c1,t) and hasnot(c1,x), c1*psolution(b,m,c,c2,y,t,n), if(hasnot(c2,t) and hasnot(c2,x), c2*psolution(b,m,c,c1,y,t,n),
	if(has(c1,t) and hasnot(c1,x), int(exp(-c*z/b)*replace(c1,t,t-x/b+z/b)*replace(c2,x,z),z,m,1,x)*exp(c*x/b)/b )));

#psolution(b_,m_,c_,ax_*t,y_,tx,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_,tx,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_,tx,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_,tx,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_,tx,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_,tx,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_,tx,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_,tx,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_,tx,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_,tx,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_,tx,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_,tx,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_,tx,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_,tx,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_,tx,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_,tx,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_,tx,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_,tx,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_,tx,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/b,y,x,m)*t-b*if(abs(c)==1,2sgn(c),sgn(c)/c/c)  )));
psolution(b_,m_,c_,t*x,y_,tx,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/c^2-(2+c*t)*b/c^3 )));
psolution(b_,m_,c_,exp(t+x),y_,tx,p_):= if(isconstant(b,c), if(1+b-c==0, exp(t+x)*x, exp(t+x)/(1+b-c) ));
psolution(b_,m_,c_,exp(t-x),y_,tx,p_):= if(isconstant(b,c), if(1-b-c==0 and p==1, exp(t-x)*x/b, exp(t-x)/(1-b-c) ));

#psolution(b_,m_,0,t*x,y_,tx,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_,tx,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_,tx,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_,tx,p_):= if(isconstant(b), if(p>1, x^2*t/2/b,int(t,t,p)*x-b*int(t,t,2p) ));

psolution(b_,m_,c_,c1_*c2_*c3_,y_,x_,n_):=if(isconstant(c1), c1*psolution(b,m,c,c2*c3,y,x,n), if(isconstant(c2), c2*psolution(b,m,c,c1*c3,y,x,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_,tx,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,y,x,m),psolution(b,m,c,c2,y,tx,n) ), 
	if(has(c1,x) and hasnot(c1,t), psolution(c/b,c1/b,y,x,m)+if(has(c2,t) and hasnot(c2,x), psolution(c,c2,y,t,n),psolution(b,m,c,c2,y,tx,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,y,x,m),psolution(b,m,c,c1,y,tx,n) ), 
	if(has(c2,x) and hasnot(c2,t), psolution(b,m,c,c1,y,tx,n)+psolution(c/b,c2/b,y,x,m), 
	psolution(b,m,c,c1,y,tx,n)+psolution(b,m,c,c2,y,tx,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_,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(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_/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(isconstant(c,d),-d/c*x^2 );
psolution(1/x_,1,x_^(-2),d_,y_,x_,2):=if(isconstant(d),-d*x^2 );
psolution(x_^(-1),1,x_^(-1),d_,y_,x_,2):=if(isconstant(d),d-d*x);
psolution(x_^m_,1,x_^m_,x_^n_,y_,x_,q_):=if(q>1 and n-m==1,1-x);

psolution(1/x_,1,x_^(-2),x_^n_,y_,x_,2):=x^(2+n)/((n+2)^2-2*(n+2)-1);
psolution(1/x_,1,x_^(-2),x_,y_,x_,2):=x^3/2;
psolution(1/x_,1,1,x_^2,y_,x_,2):= -x^2;
psolution(1/x_,1,1/x_,x_^n_,y_,x_,2):=	if(n>0,-sum((1+n)!/(1+n-k)!*x^(n+1-k),k,0,n-1) );
#psolution(1/x,1,1/x,x^2,y,x,2):= -3x^2-x^3;
#psolution(1/x,1,1/x,x^3,y,x,2):= -24x^2-8*x^3-x^4;
#psolution(1/x,1,1/x,x^4,y,x,2):= -360*x^2-120*x^3-15*x^4-x^5;

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(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(x_,1,-x_,exp(x_),y_,x_,2):=exp(x);
psolution(-x_,1,x_,exp(x_),y_,x_,2):=exp(x);

psolution(b_,1,0,x_^2,y_,x_,n_):= if(isconstant(b) and n>2, -x/b-x^3/b/3 );
psolution(b_,1,0,x_,y_,x_,n_):= if(isconstant(b) and n>2, -1/b/2*x^2 );

psolution(b_,1,0,x_^2,y_,x_,2):= if(isconstant(b), -2x/b-x^2/b-x^3/b/3 );
psolution(b_,1,0,x_,y_,x_,2):= if(isconstant(b), -1/b/2*x^2-x/b^2 );
psolution(b_,1,0,exp(b_*x_),y_,x_,2):= exp(b*x)*x/b;
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(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);
psolution(b_,mm_,c_,c_*x_^m_,y_,x_,n_):=if(isinteger(m) or isinteger(n-m) and n>m and m>0 and mm>1, -x^m);
psolution(b_,m_,c_,c_*x_,y_,x_,n_):= if(m>1, -x);
psolution(b_,m_,c_,c_,y_,x_,n_):= if(m>0, -1);


#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) 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)) ));
psolution2(c_,d_,y_,x_,n_):=if(isfree(d/c,x), -d/c,  if(n<=1,block(dsolve:=n!*integrate(c,x,n), integrate(d/mittag(n,dsolve),x,n)*mittag(n,dsolve)) ));
#psolution2(b_,x_^n_, y_,x_,q_) := If(isfree(b,x) and isfree(b,y) and n>0, if(q==n,-x^n/b-n!, if(q-n>0,-x^n/b, if(isinteger(n) or isinteger(n-q),-sum(b^(-1-k/q)*n!/(n-k)!*x^(n-k),k,0,n,q),-gamma(n+1,b*x)*exp(b*x) ))));

psolution2(c_,d_,y_,x_,2):=if(isfree(c,x),if(isfree(d,x), -d/c,  
	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))/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)) )));
psolution2(c_*x_,f_,y_,x_,2):= if(isconstant(c),pi*Bi(cbrt(c)*x)*int(Ai(cbrt(c)*x)*f,x)/cbrt(c)-pi*Ai(cbrt(c)*x)*int(Bi(cbrt(c)*x)*f,x)/cbrt(c));
psolution2(x_,f_,y_,x_,2):= pi*Bi(x)*int(Ai(x)*f,x)-pi*Ai(x)*int(Bi(x)*f,x);
psolution2(x_^(-1),d_,y_,x_,2):= 2sqrt(x)*besselI(1,2sqrt(x))*int(besselK(1,2sqrt(x))*sqrt(x)*d,x)-2sqrt(x)*besselK(1,2sqrt(x))*int(besselI(1,2sqrt(x))*sqrt(x)*d,x);
psolution2(x_^(-2),d_,y_,x_,2):= x^(1/2 + sqrt(5)/2)* int(x^(1/2 - sqrt(5)/2) *d,x)/sqrt(5) - x^(1/2 - sqrt(5)/2)* int(x^(1/2 + sqrt(5)/2)*d,x)/sqrt(5);
psolution2(b_*x_^(-2),d_,y_,x_,2):=  if(isconstant(b),-i*x^(1/2 + 1/2 i *sqrt(-4 b - 1))* int(( x^(1/2 - 1/2 i* sqrt(-4 b - 1))* f(x)),x)/sqrt(-4 b - 1) + i*x^(1/2 - 1/2 i* sqrt(-4 b - 1))* int( (x^(1/2 i* sqrt(-4 b - 1) + 1/2)* f(x)),x)/sqrt(-4 b - 1) );
psolution2(exp(x_),d_,y_,x_,2):= 2 besselI(0,2 sqrt(exp(x)))* int(besselK(0,2 sqrt(exp(x)))* d,x) -2 besselK(0,2 sqrt(exp(x))) *int(besselI(0,2 sqrt(exp(x)))* d,x);

psolution(c_,E(n_,x_^n_),y_,x_,n_):=if(isfree(c,y) and isfree(c,x), if(c==1, 2/sqrt(pi)*E(n,x^n)*x^n, E(n,x^n)/(1-c) ));
psolution(c_,E(n_,b_*x_^n_),y_,x_,n_):=if(isfree(b,c,y) and isfree(b,c,x), if(c==b, 2/sqrt(pi)*E(n,b*x^n)*x^n,E(n,b*x^n)/(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,n,x) and isfree(b,y) and n>0, if(q==n,-x^n/b-n!, if(q-n>0,-x^n/b,  if(isinteger(n) or isinteger(n-q),-sum(b^(-1-k/q)*n!/(n-k)!*x^(n-k),k,0,n,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(b_*x_^m_,c_,y_,x_,q_) := if(isconstant(b,c),if(m+q==0,if(b==q!,c*log(x)*x^q/(fallingfactorial(4,q-1)-1),c/x^m/(q!-b)), if(m<0 and m+q>0,-c/b/x^m )));
psolution(b_/x_,c_,y_,x_,n_):= if(n>1 and isconstant(b,c),-c/b*x );
psolution(x_^m_,c_,y_,x_,q_) := if(isconstant(c),if(m+q==0,c/x^m/(q!-1), if(m<0 and m+q>0,-c/x^m )));
psolution(1/x_,c_,y_,x_,n_):= if(n>1 and isconstant(c),-x );

#psolution(b_*x_^m_,x_^n_,y_,x_,q_) := if(isconstant(b),if(m== -q,x^(n-m)/(fallingfactorial(n-m,q)-b), -1/b*x^(n-m)-if(q>=n-m,0, fallingfactorial(n-m,q)/b^q ) ));
psolution(b_*x_^m_,x_^n_,y_,x_,q_) := if(isconstant(b),if(q-n+m>0 and n-m>0,-x^(n-m)/b, if(q==n-2m,-x^(n-m)/b-fallingfactorial(n-m,q)/b^q, if(m== -q,x^(n-m)/(fallingfactorial(n-m,q)-b)  ))));
psolution(b_*x_^m_,x_,y_,x_,q_) := if(m<0 and q+m>=1,-x^(1-m)/b, if(m== -q,x^(1+q)/((1+q)!-b)  ));
#psolution(c_*x_,x_^n_,y_,x_,q_):= if(isconstant(c), if(q>=n-1,-1/c*x^(n-1), if(n-q==2,-(n-1)!/c/c-x^3/c, if(q==2,pi*Bi(cbrt(c)*x)*int(Ai(cbrt(c)*x)*x^n,x)/cbrt(c)-pi*Ai(cbrt(c)*x)*int(Bi(cbrt(c)*x)*x^n,x)/cbrt(c) ))));
#psolution(c_*x_,x_^n_,y_,x_,q_):= if(isconstant(c), if(q-n>= -1,-1/c*x^(n-1), if(n-q==2,-(n-1)!/c^2-x^3/c )));
psolution(x_^m_,x_^n_,y_,x_,q_) := if(q-n+m>0 and n-m>0,-x^(n-m), if(q==n-2m,-x^(n-m)-fallingfactorial(n-m,q), if(m+q==0,x^(n-m)/(fallingfactorial(n-m,q)-1), if(q>1 and n>0,-sum(n!/(n-k)!*(n+1)!/(n+1-k)!*(n+2)!/(n+2-k)!*x^(n-m-k),k,0,n,q+m)  ))));
psolution(x_^m_,x_,y_,x_,q_) := if(m<1 and q+m>=1,-x^(1-m), if(m== -q,x^(1-m)/((1-m)!-1) ));
#psolution(x_,x_^n_,y_,x_,q_):= if(q-n== -2,-(n-1)!-x^(n-1), if(q-n== -1,-(n)!*x-x^(n-1), if(q-n>-1,-x^(n-1), if(q==2,pi*Bi(x)*int(Ai(x)*x^n,x)-pi*Ai(x)*int(Bi(x)*x^n,x) ))));
psolution(b_*x_,x_^n_,y_,x_,q_):= if(isconstant(b), if(n>1 and n-q<1,-x^(n-1)/b ,if(q>1 and n-q>1 and  n<2q+2,-fallingfactorial(n-1,q)*x^(n-2-q)/b^2-x^(n-1)/b )));
#psolution(x_,x_^n_,y_,x_,q_):= if(n>1 and n<1+q,-x^(n-1), if(q>1 and n-q>1 and n<2q+2,-fallingfactorial(n-1,q)*x^(n-2-q)-x^(n-1),  if(q>1 and n>1,-sum((n-1)!*(n)!/(n-1-k)!/(n-k)!*x^(n-1-k),k,0,n-1,q+1) )));
#psolution(1/x_,x_^n_,y_,x_,q_):= if(q==n,-(1+n)!*x^2-x^(1+n), if(q-n==1,-q!*x-x^q, if(n> -1 and q-n>1,-x^(1+n), if(q>1 and n>0,-sum((n+1)!/(n-2/3*k+1)!*x^(n-k+1),k,0,n+1,q-1) ))));
psolution(b_/x_,x_^n_,y_,x_,q_):= if(isconstant(b),if(q>2 and q-n>=0 and q-n<=1,-fallingfactorial(n+1,q)*x^(n+2-q)/b^2-x^(1+n)/b, if(n> -1 and q-n>1,-x^(1+n)/b, if(q>1 and n>0,-sum(n!*(n+1)!/(n+1-k)!/(n-k)!/b^(1+k)*x^(n+1-k),k,0,n,q-1) ))));
#psolution(1/x_,x_^n_,y_,x_,q_):= if(q>2 and q-n>=0 and q-n<=1,-fallingfactorial(n+1,q)*x^(n+2-q)-x^(1+n), if(n> -1 and q-n>1,-x^(1+n), if(q==2 and n>0,-sum((n)!/(n-k)!*(n+1)!/(n+1-k)!*x^(n+1-k),k,0,n,q-1) )));
#psolution(1/x_,x_^n_,y_,x_,q_):= if(q>2 and n/(q-1)<3 and q-n>= -2 and q-n<=1,-fallingfactorial(n+1,q)*fallingfactorial(n+2-q,q)*x^(n+3-2q)-fallingfactorial(n+1,q)*x^(n+2-q)-x^(1+n), ;
psolution(1/x_,x_^n_,y_,x_,q_):= if(n> -1,
	if(n<q-1,-x^(1+n), 
	if(n<2*q-2,-fallingfactorial(n+1,q)*x^(n+2-q)-x^(1+n), 
	if(n<3*q-3,-fallingfactorial(n+1,q)*fallingfactorial(n+2-q,q)*x^(n+3-2q)-fallingfactorial(n+1,q)*x^(n+2-q)-x^(1+n)
	))));
psolution(x_,x_^n_,y_,x_,q_):= if(n>1,
	if(n<q+1,-x^(n-1), 
	if(n<2*q+2,-fallingfactorial(n-1,q)*x^(n-2-q)-x^(n-1), 
	if(n<3*q+3,-fallingfactorial(n-1,q)*fallingfactorial(n-2-q,q)*x^(n-3-2q)-fallingfactorial(n-1,q)*x^(n-2-q)-x^(n-1)
	))));
	
psolution(1/x_,x_^n_,y_,x_,2):=	if(n>0,-sum(n!/(n-k)!*(n+1)!/(n+1-k)!*x^(n+1-k),k,0,n) );
psolution(2*x_^(-2),c_,y_,x_,2) := if(isconstant(c),c*log(x)*x^2/3 );
psolution(6*x_^(-3),c_,y_,x_,3) := if(isconstant(c),c*log(x)*x^3/11 );
psolution(24*x_^(-4),c_,y_,x_,4):= if(isconstant(c),c*log(x)*x^4/50 );
psolution(c_,sin(x_),y_,x_,2):= if(c== -1,-cos(x)*x/2, if(isfree(c,x), -sin(x)/(c + 1) ));
psolution(c_,cos(x_),y_,x_,2):= if(c== -1,sin(x)*x/2, if(isfree(c,x), -cos(x)/(c + 1) ));
psolution(c_,tan(x_),y_,x_,2):= if(c== -1,-cos(x)*log(sec(x)+tan(x)) );
psolution(c_,cot(x_),y_,x_,2):= 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_,sin(x_)*x_,y_,x_,2):= if(c== -1,-cos(x)*x^2/4+sin(x)*x/4, if(isfree(c,x), c/2*cos(x)/(c+1)^2-sin(x)*x/(c + 1) ));
psolution(c_,cos(x_)*x_,y_,x_,2):= if(c== -1,sin(x)*x^2/4+cos(x)*x/4, if(isfree(c,x), -c/2*sin(x)/(c+1)^2-cos(x)*x/(c + 1) ));
#psolution(x_,x_^7,y_,x_,2):= -180-30*x^3-x^6;
psolution(1/x_,x_,y_,x_,2):= -2x-x^2;
psolution(1/x_,x_,y_,x_,1.5):= -3x+(-2.256758334191025)*x^1.5-x^2;
psolution(1/x_,x_^1.5,y_,x_,1.5):=(-9.970052911343528)*x+(-7.5)*x^1.5-3.3233509704478426*x^2-x^2.5;
psolution(1/x_,x_^6,y_,x_,3):= (-75600.0)*x-12600*x^3-210*x^5-x^7;

#psolution(1/x_,x_^4,y_,x_,3):= -360x-60x^3-x^5;
#psolution(1/x_,x_^4,y_,x_,2):= (-2880.0)*x-1440x^2-240*x^3-20*x^4-x^5;
#psolution(1/x_,x_^3,y_,x_,2):= (-144.0)*x+(-72.0)*x^2-12*x^3-x^4;
#psolution(1/x_,x_^2,y_,x_,2):= -12x-6x^2-x^3;
#psolution(1/x_,x_^3.5,y_,x_,2.5):=(-177.1875)*x^1.5-26.171388892276763*x^3-x^4.5;
#psolution(1/x_,x_^5,y_,x_,2.5):=(-10967.845504168383)*x^1.5-1620.0*x^3.0-61.8996571663824*x^4.5-x^6;
#psolution(1/x_,x_^4.5,y_,x_,2.5):=(-2878.8527781504436)*x-866.25*x^2.5-47.98087963584073*x^4.0-x^5.5;

#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(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) ));

#psolution2(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(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;
psolution(b_,c_*d_,y_,t,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_,y_,x_,n_):=psolution(c,c1,y,x,n)+psolution(c,c2,y,x,n);

psolution(c_,b_*c_,y_,x_,n_):=if(isfree(b,x),  -b);
psolution(a_*c_,c_,y_,x_,n_):=if(isfree(a,x),  -1/a);
psolution(a_*c_,b_*c_,y_,x_,n_):=if(isfree(a,b,x), -b/a);
psolution(c_,c_*x_^m_,y_,x_,n_):=if(isinteger(m) or isinteger(n-m) and n>m and m>0, -x^m);
psolution(c_,c_*x_,y_,x_,n_):=if(n>1, -x);
psolution(a_*x_,a_*x_^n_,y_,x_,p_):=if(p>n-1 and isinteger(n),-x^(n-1));
psolution(a_*x_^m_,a_*x_^n_,y_,x_,p_):=if(p>n-m and isinteger(n-m),-x^(n-m));
psolution(c_,c_,y_,x_,n_):= -1;

psolution(0,c_,y_,x_,n_):=if(isfree(c,y), int(c,x,n) );