#	differential equation solver;
#	solve differential equation and fractional differential equation;
#	e.g. dsolve( y' = x*y, y,x);
#	dsolve( ds(y,x,0.5)=y^2 );

#dsolve(a_,y_,x_) := If(isfree(a,y),integrate(a,x), If(isfree(a/y,y),exp(integrate(a/y,x)), If(isfree(a,x), (integrate(1/a, y)=x) )));

dsolve(y_,y_,x_) := exp(x);
dsolve(f(a_*x_+b_*y_),y_,x_) := -a/b*x+inversef(-a/b)/b;
dsolve(f(x_+b_*y_),y_,x_) := -1/b*x+inversef(-1/b)/b;
dsolve(f(a_*x_+y_),y_,x_) := -a*x+inversef(-a);
dsolve(f(x_+y_),y_,x_) := -x+inversef(-1);
dsolve(sin(a_*x_+b_*y_),y_,x_) := -a/b*x+asin(-a/b)/b;
dsolve(cos(a_*x_+b_*y_),y_,x_) := -a/b*x+acos(-a/b)/b;
dsolve(sin(a_*x_+b_*y_+C_),y_,x_) := If(isfree(a,b,c,x) and isfree(a,b,c,y),-a/b*x+asin(-a/b)/b-c/b);
dsolve(cos(a_*x_+b_*y_+c_),y_,x_) := If(isfree(a,b,c,x) and isfree(a,b,c,y),-a/b*x+acos(-a/b)/b-c/b);

dsolve(y'=a_,y,x_) :=dsolve(a,y,x);
dsolve(y'=y,y,x_) := exp(x);

dsolve(a_+b_, y_,x_) := If(isfree(a,b,y), integrate(a+b,x), 
	If(has(a,y) and has(b,y) and isfree(d(a+b,y),y), exp(integrate(d(a+b,y),x)), 
	If(isfree(b,y) and isfree(d(a,y),y),  	block(dsolve:=int(d(a,y),x),exp(dsolve)*int(b*exp(-dsolve),x)+exp(dsolve)), 
	If(isfree(a,y) and isfree(d(b,y),y), 	block(dsolve:=int(d(b,y),x),exp(dsolve)*int(a*exp(-dsolve),x)+exp(dsolve)) ))));
dsolve(a_*y_+b_*y_^n_,y_,x_):=If(isfree(a,b,n,y) , block(dsolve:=integrate(a,x), exp(dsolve)*(1+(1-n)*integrate(b*exp((n-1)*dsolve),x))^(1/(1-n)) ) );
dsolve(a_*y_+y_^n_*z_,y_,x_):=If(isfree(a,z,n,y) , block(dsolve:=integrate(a,x), exp(dsolve)*(1+(1-n)*integrate(z*exp((n-1)*dsolve),x))^(1/(1-n)) ) );
dsolve(y_+y_^n_,y_,x_):=If(isfree(n,y),  exp(x)*(1-exp((n-1)*x))^(1/(1-n)) );

#dsolve(a_+b_*y_,y_,x_):=If(isfree(a,b,y), block(dsolve:=integrate(b,x), (exp(dsolve)*integrate(a*exp(-dsolve),x)+exp(dsolve)) ) );
#dsolve(a_+y_*z_,y_,x_):=If(isfree(a,z,y), block(dsolve:=integrate(z,x), (exp(dsolve)*integrate(a*exp(-dsolve),x)+exp(dsolve)) ) );
#dsolve(a_+b_*y_,y_,x_):=If(isfree(a,b,x) and isfree(a,b,y), exp(b*x)-a/b);
#dsolve(a_+y_,y_,x_) := If(isfree(a,y),  -integrate(a,x), If(isfree(a,x), solve(integrate(1/a, y)= -x,y) ));
#dsolve(y_'+z_,y_,x_) := If(isfree(z,y),  -integrate(z,x), If(isfree(z,x), solve(integrate(1/z, y)= -x,y) ));
#dsolve(z_+y_, y_,x_) := If(isfree(z,y), exp(x)-z);

dsolve(z_+y_*z_,y_,x_):=If(isfree(z, y),y=exp(integrate(z,x))-1);
dsolve(a_+a_*y_,y_,x_):=If(isfree(a, y),y=exp(integrate(a,x))-1);

dsolve(x_+y_, y_,x_) := -x-1;
dsolve(a_*x_+y_, y_,x_) := If(isfree(a, x) and isfree(a, y), -a*x-a);
dsolve(x_+b_*y_,y_,x_)  := If(isfree(b, x) and isfree(b, y), -1/b*x-1/b^2);
dsolve(a_*x_+b_*y_,y_,x_):=If(isfree(a,b, x) and isfree(a,b, y), -a/b*x-a/b^2);
dsolve(x_+y_+c_, y_,x_) := If(isfree(c,x) and isfree(c,y), -x-1-c);
dsolve(x_+b_*y_+c_, y_,x_) := If(isfree(b,c,x) and isfree(b,c,y), -1/b*x-1/b^2-c/b);
dsolve(a_*x_+y_+c_, y_,x_) := If(isfree(a,c,x) and isfree(a,c,y), -a*x-a-c);
dsolve(a_*x_+b_*y_+c_, y_,x_) := If(isfree(a,b,c,x) and isfree(a,b,c,y), -a/b*x-a/b^2-c/b);

dsolve(g_*y_, y_,x_) := If(isfree(g,y), exp(integrate(g,x)));
dsolve(y_*z_, y_,x_) := If(isfree(z,y), exp(integrate(z,x)));
dsolve(g_*y_^n_, y_,x_) := If(isfree(g,y), ((1-n)*integrate(g,x))^(1/(1-n)) );

dsolve(ds(y_,x_,q_)=a_,y_,x_) := If(isfree(a,y),d(a,x,-q), If(isfree(a/y,y),mittag(q,d(a/y,x,-q)), If(isfree(a,x), solve(d(1/a, y,-q)=x^q/q!,y) )));
dsolve(ds(y_,x_,q_)=a_+b_, y_,x_) := If(isfree(a,b,y), d(a+b,x,-q), 
	If(isfree(a,b,x), solve(d(1/(a+b),y,-q)=x^q/q!,y), 
	If(isfree(a,y) and isfree(d(b,y),y), 	block(dsolve:=d(d(b,y),x,-q),simplify(mittag(q,dsolve)*d(a*mittag(q,-dsolve),x,-q))+mittag(q,dsolve)), 
	If(isfree(b,y) and isfree(d(a,y),y), 	block(dsolve:=d(d(a,y),x,-q),simplify(mittag(q,dsolve)*d(b*mittag(q,-dsolve),x,-q))+mittag(q,dsolve)) ))));

dsolve(ds(y_,x_, q_)=a_+y_, y_,x_) := If(isfree(a,y), If(isfree(a,x), If(q>0,exp(x)-a, a*mittag(-q,x^(-q)) ), simplify(mittag(q,x^q/q!)*d(a*mittag(q,-x^q/q!),x,-q))+mittag(q,x^q/q!) ));
dsolve(ds(y_,x_, q_)=a_+b_*y_, y_,x_) := If(isfree(a,b,x) and isfree(a,b,y), If(q<0, a*mittag(-q,1/b*x^(-q)), If(b>0, exp(b^(1/q)*x)-a/b, mittag(q,b*x^q/q!)-a/b  )), mittag(q,d(b,x,-q))-a);

#dsolve(ds(y_,x_, q_)=x_^q_+y_, y_,x_) := -x^q-q!;
#dsolve(ds(y_,x_, q_)=x_^n_+b_*y_, y_,x_) := 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)) )));

#dsolve(ds(y_,x_, q_)=x_^n_+b_*y_, y_,x_) := If(isfree(b,x) and isfree(b,y), If(n==q, -x^n/b-q!/b^2, n!*x^(n+q)*mittag(q,1+n+q,b*x^q) ));
#dsolve(ds(y_,x_, q_)=a_*x_^n_+b_*y_, y_,x_) := If(isfree(a,b,x) and isfree(a,b,y), If(n==q, -a*x^n/b-a*q!/b^2, a*n!*x^(n+q)*mittag(q,1+n+q,b*x^q) ));
dsolve(ds(y_,x_, q_)=x_^n_+y_, y_,x_) := If(n==q, -x^n-n!, n!*x^(n+q)*mittag(q,1+n+q,x^q));
dsolve(ds(y_,x_, q_)=a_*x_^n_+y_, y_,x_) := If(isfree(a,x) and isfree(a,y), If(n==q, -a*x^n-a*n!, a*n!*x^(n+q)*mittag(q,1+n+q,x^q) ));

dsolve(ds(y_,x_, q_)=x_^n_+y_+c_, y_,x_) := If(isfree(c,x) and isfree(c,y), If(n==q, -x^n-n!-c, n!*x^(n+q)*mittag(q,1+n+q,x^q)-c));
dsolve(ds(y_,x_, q_)=a_*x_^n_+y_+c_, y_,x_) := If(isfree(a,c,x) and isfree(a,c,y), If(n==q, -a*x^n-a*n!-c, a*n!*x^(n+q)*mittag(q,1+n+q,x^q)-c ));
#dsolve(ds(y_,x_, q_)=x_^n_+b_*y_+c_, y_,x_) := If(isfree(b,c,x) and isfree(b,c,y), If(n==q, -x^n/b-q!/b^2-c/b, n!*x^(n+q)*mittag(q,1+n+q,b*x^q)-c/b ));
#dsolve(ds(y_,x_, q_)=a_*x_^n_+b_*y_+c_, y_,x_) := If(isfree(a,b,c,x) and isfree(a,b,c,y), If(n==q, -a*x^n/b-a*q!/b^2-c/b, a*n!*x^(n+q)*mittag(q,1+n+q,b*x^q)-c/b ));

dsolve(ds(y_,x_, q_)=a_*x_+b_*y_, y_,x_) := If(isfree(a,b,x) and isfree(a,b,y), a*x^(1+q)*mittag(q,2+q,b*x^q));
dsolve(ds(y_,x_, q_)=a_*x_+y_, y_,x_) := If(isfree(a,x) and isfree(a,y), a*x^(1+q)*mittag(q,2+q,x^q));
dsolve(ds(y_,x_, q_)=x_+b_*y_, y_,x_) := If(isfree(b,x) and isfree(b,y), x^(1+q)*mittag(q,2+q,b*x^q) );
dsolve(ds(y_,x_, q_)=x_+y_, y_,x_) :=  x^(1+q)*mittag(q,2+q,x^q);
dsolve(ds(y_,x_, q_)=z_+y_*z_,y_,x_):=If(isfree(z, y),mittag(q,d(z,x,-q))-1);
dsolve(ds(y_,x_, q_)=a_+a_*y_,y_,x_):=If(isfree(a, y),mittag(q,d(a,x,-q))-1);
dsolve(ds(y_,x_, q_)=a_*x_+y_+c_, y_,x_) := If(isfree(a,c,x) and isfree(a,c,y), If(q>0,a*x^(1+q)*mittag(q,2+q,x^q)-c,a*x^(1+q)*mittag(q,2+q,x^q) ));
dsolve(ds(y_,x_, q_)=x_+b_*y_+c_, y_,x_) := If(isfree(b,c,x) and isfree(c,b,y), If(q>0,x^(1+q)*mittag(q,2+q,b*x^q)-c/b,x^(1+q)*mittag(q,2+q,b*x^q) ));
dsolve(ds(y_,x_, q_)=x_+y_+c_, y_,x_) := If(isfree(c,x) and isfree(c,y), If(q>0,x^(1+q)*mittag(q,2+q,x^q)-c,x^(1+q)*mittag(q,2+q,x^q) ));
dsolve(ds(y_,x_, q_)=a_*x_+b_*y_+c_, y_,x_) := If(isfree(a,b,c,x) and isfree(a,b,c,y), If(q>0, a*x^(1+q)*mittag(q,2+q,b*x^q)-c/b,a*x^(1+q)*mittag(q,2+q,b*x^q) ));
dsolve(ds(y_,x_, q_)=a_*x_+x_*y_, y_,x_) := If(isfree(a,x) and isfree(a,y), mittag(q,x^(1+q)/(1+q)!)-a);
dsolve(ds(y_,x_, q_)=a_*x_+b_*x_*y_, y_,x_) := If(isfree(a,b,x) and isfree(a,b,y), mittag(q,b*x^(1+q)/(1+q)!)/b-a/b);

dsolve(ds(y_,x_, q_)=b_*y_^n_, y_,x_) := If(isfree(b,y), If(isfree(b,x), b^(1-n)*fallingfactorial(q/(1-n),q)^(1/(n-1))*x^(q/(1-n)), (1/fallingfactorial(-n,q)*d(b,x,-q))^(1/(q-n)) ));
dsolve(ds(y_,x_, q_)=y_^n_, y_,x_) := fallingfactorial(q/(1-n),q)^(1/(n-1))*x^(q/(1-n)) ;

#dsolve(ds(y_,x_, q_)=x_*y_+a_, y_,x_) :=  exp(q/(1+q)*x^(1/q+1))-a;

dsolve(ds(y_,x_, q_)=b_*y_, y_,x_) := If(b>0, exp(integrate(b^(1/q),x)), If(b<0, mittag(q,b*x^q/q!), mittag(q,d(b,x,-q)) ));
dsolve(ds(y_,x_, q_)=b_*x_*y_, y_,x_) := If(b>0, exp(b^(1/q)*q/(1+q)*x^(1/q+1)), If(b<0, If(q>0, mittag(q,b*x^(q+1)/(1+q)!),mittag(-q,1/b*x^(1-q)) )));
dsolve(ds(y_,x_, q_)=y_*z_, y_,x_) := If(isfree(z,y), mittag(q,d(z,x,-q)) );
dsolve(ds(y_,x_, q_)=x_*y_, y_,x_) :=  If(q>0,exp(q/(1+q)*x^(1/q+1)), mittag(-q,x^(1-q)/(1-q)!));
dsolve(ds(y_,x_, q_)=y_, y_,x_) := exp(x);

dsolve(f_,y) := dsolve(y'-f,y,x);
dsolve(f_,y(x_)) := dsolve(replace(d(y(x),x)-f,y(x),y),y,x);
dsolve(a_=b_,y_):=dsolve(a-b,y);
dsolve(y'=f_,y_) := dsolve(f,y,x);
dsolve(a_+y',y_) := dsolve(-a,y,x);
dsolve(y'+z_,y_) := dsolve(-z,y,x);
dsolve(a_+y'+z_,y_) := dsolve(-a-z,y,x);
dsolve(y'+z_+zz_,y_) := dsolve(-zz-z,y,x);
dsolve(d(y(x_),x_)=f_,y(x_)) := dsolve(replace(f,y(x),y),y,x);
dsolve(ds(y_)=a_,y_) := dsolve(a,y,x);
dsolve(ds(y_,x_)=a_,y_) := dsolve(a,y,x);
dsolve(ds(y_,x_,q_)=a_,y_):= dsolve(ds(y,x,q)=a,y,x);
dsolve(a_+ds(y_,x_,q_),y_):= dsolve(ds(y,x,q)= -a,y,x);

dsolve(a_):=dsolve(a,y);