@Joe Riel 

Yes, you are right.

@Mikhak 

Here symbolic solutions does not make much sense. Try:

solve(F, x, allsolutions);
F1:=convert(F, rational);
solve(F1, x, allsolutions);

It seems that you want numeric solutions.

f := (1 + a * cos(theta1*x)) * cos(theta2*x) - b:
F:=eval(f, [theta1=2, theta2=24.5, a=0.25 , b=0.1]):
RootFinding:-Analytic( F, x, re=0..4*Pi, im=-0.1 .. 0.1); sort([%]);

[0.0608404851442370, 0.195661496064193, 0.317169250177469, 0.452338844731748, 0.573321271654020, 0.709193982762240, 0.829303334433960, 0.966198716646160, 1.08517751285740, 1.22323304796064, 1.34112176470098, 1.48005365921274, 1.59739755110209, 1.73644235409529, 1.85412330733945, 1.99244553947813, 2.11113047152258, 2.24831526005946, 2.36815894926479, 2.50426092567311, 2.62505648897743, 2.76036869003974, 2.88177847976375, 3.01665365639252, 3.13832384040376, 3.27311241579994, 3.39469611285126, 3.52974541499388, 3.65089259774324, 3.78655671358430, 3.90691457307736, 4.04353049236051, 4.16280435532856, 4.30057605523926, 4.41870733647864, 4.55747629551925, 4.67487862383864, 4.81398122172356, 4.93149765193513, 5.07006063884640, 5.18846029384880, 5.32593928588975, 5.44550314841780, 5.58185358677780, 5.70244066316500, 5.83791777406930, 5.95920741126850, 6.09415862817895, 6.21579641986160, 6.35057419449755, 6.47221198618020, 6.60716320309065, 6.72845284028987, 6.86392995119415, 6.98451702758135, 7.12086746594140, 7.24043132846945, 7.37791032051040, 7.49630997551275, 7.63487296242400, 7.75238939263560, 7.89149199052055, 8.00889431883995, 8.14766327788050, 8.26579455911990, 8.40356625903060, 8.52284012199870, 8.65945604128180, 8.77981390077485, 8.91547801661595, 9.03662519936530, 9.17167450150790, 9.29325819855925, 9.42804677395540, 9.54971695796665, 9.68459213459545, 9.80600192431945, 9.94131412538175, 10.0621096886861, 10.1982116650944, 10.3180553542997, 10.4552401428366, 10.5739250748805, 10.7122473070197, 10.8299282602638, 10.9689730632571, 11.0863169551464, 11.2252488496582, 11.3431375663985, 11.4811931015018, 11.6001718977130, 11.7370672799252, 11.8571766315982, 11.9930493427051, 12.1140317696274, 12.2492013641817, 12.3707091182950, 12.5055301292148]

You have all the 98 real solutions in the interval 0 .. 4*Pi.
F is 4*Pi - periodic, so any other real solution can be obtained by adding an integer multiple of 4*Pi.

@mmcdara 

A Maple command to squeeze polynomials is compoly. Unfortunately it fails for almost any random polynomial.
My example was not presented as a challenge (and the simplified form was given).
Here is also a simple one:

F := 1238*a^9 + (6927*b + 3756)*a^8 + (40476*b^2 + 26226*b + 5334)*a^7 + (145584*b^3 + 144582*b^2 + 46947*b + 5376)*a^6 + (326904*b^4 + 436890*b^3 + 218880*b^2 + 48384*b + 4032)*a^5 + (489558*b^5 + 815670*b^4 + 543885*b^3 + 181440*b^2 + 30240*b + 2016)*a^4 + (489988*b^6 + 979926*b^5 + 816510*b^4 + 362880*b^3 + 90720*b^2 + 12096*b + 672)*a^3 + (314936*b^7 + 734882*b^6 + 734877*b^5 + 408240*b^4 + 136080*b^3 + 27216*b^2 + 3024*b + 143)*a^2 + (118084*b^8 + 314894*b^7 + 367396*b^6 + 244944*b^5 + 102060*b^4 + 27216*b^3 + 4536*b^2 + 430*b + 16)*a + 19686*b^9 + 59055*b^8 + 78735*b^7 + 61236*b^6 + 30618*b^5 + 10206*b^4 + 2268*b^3 + 323*b^2 + 25*b;

It can be written (in 1D) using 534 characters (using simplify/size), i.e.  length(convert(F,string)) = 534.
There is a simplified form using only 49 characters.
It would be very nice to have in Maple such a "squeezer".

I did not remove it! Probably somebody else did.

I just asked why it was transformed into a question, from its initial place.

In all cases, it would be nice that any such action be mentioned by the person who decides to do so.

I wonder why this comment @chandrashekhar was converted into a question.

@Rouben Rostamian  

Very nice, thank you.
So, it seems that pdsolve accepts discontinuous bc's but rejects discontinuous pdes, and converting the 2nd order PDE into a 1st order system does not always work.
The solution is to approximate the discontinuous pde or to implement a custom finite difference method, as you did.
[There should be some pdsolve options, but they seem to be hard to find].

The pde approximation method (also yours!) works well for your last example:

restart;
uexact := piecewise(x<0, erf(x/(2*sqrt(t))), 1/2*erf(x/(4*sqrt(t)))):
#K := x -> piecewise(x<=0, 1,4):
K := x -> 1+3*(1+erf(100*x))/2:
pde := { diff(u(x,t),t) =  diff(K(x)*diff(u(x,t),x), x)}:
bc:={ u(-1,t)=-erf(1/(2*sqrt(t))), u(1,t)=erf(1/(4*sqrt(t)))/2, u(x,0)=piecewise(x<0, -1, +1/2)}:
sol:=pdsolve(pde, bc, numeric, spacestep=1/200, timestep=0.0001):
p1:=sol:-animate(u, t=0..0.1, frames=25, color=red):
p2:=plots:-animate(plot,[uexact, x=-1..1, color=blue], t=0..0.1, frames=25):
plots:-display(p1,p2);

@Rouben Rostamian

It would be nice to know an exact solution for an example where c(x) (=K(x)) is also discontinuous.
It is in this case where I don't know theoretical results.

@Bart 

Yes, it seems that pdsolve rejects not only a discontinuous u(x,0) but also a continuous but non-smooth one.
Rouben's solution accepts any u(x,0) because it does not check it. It just solves the generated linear system.
I do not know what the theory says for such conditions and about the nature of the solution (when it exists).

(for Example 2)

restart;
K := x -> piecewise(x<0.5, 1, 10):
pde := {diff(u(x,t),t)=diff(v(x,t),x),
        v(x,t)=K(x)*diff(u(x,t),x)}: 
bc:={ u(0,t)=0, u(1,t)=1, u(x,0)=2 }:
sol:=pdsolve(pde, bc, numeric, spacestep=1/200):
sol:-animate(u, t=0..0.15);

@federicie

f:=proc(a,b,ev)
   local i; global A,B,F;
   if ev<>true then return 'last'(A,B)=F fi;
   A,B:=a,b;
   F:=add(evalf(sin(a+b+i)), i=1..10^5)
end:

@Carl Love 
In the 2019 help page:

This feature is currently not supported in 2-D math input

@Carl Love 

Maybe I am too conservative but I prefer the old style

ur:= proc(r) local i:=op(procname); A[i]*r + B[i]/r end;

@WebSuccess 

Unfortunately you have not checked your statement.

@byrktr 

Sorry, I don't know. But usually in such problems some heuristic is used.
BTW, I have corrected a typo in cond3.

Edit.
Here is an idea.
Solve the problem in real x[i,j,k],  in the interval 0..1. Maple can do it quickly.
One obtains an "approx" solution (with some conditions violated). Maybe you can start from here.

i1:=indets([cond1])<=~1:
i0:=indets([cond1])>=~0:
i01:=i0[],i1[]:

fracsol:=Optimization:-LPSolve( add(q[i],i=1..20), {seq(cond||i, i=1..6), i01} );
sol01_approx:=map(t -> (lhs(t)=round(rhs(t))), fracsol[2]);