@sand15 The solution space for msolve(..., n) is the set of equivalence classes of remainders mod n, which is a finite set with n elements. When working with modular arithmetic at anything other than an elementary level, it goes without saying that objects that are written as ordinary integers are actually representatives of their equivalence classes. In other words, x means {x + k*n | k in Z}. This a totally standard notation in modular arithmetic. 

IMO, algsubs is junk, and it should be deprecated. I would only use it as a last resort. Indeed, I haven't used it in at least 16 years. At the very least, its help page should advise that there are several newer and better commands, such as simplify with side relations. 

As you said, you're using an approximation method, so of course it deviates from the exact solution. In particular, approximation methods for IVPs tend to deviate more the further you get from the initial point.

So, what kind of help are you expecting? There's nothing wrong with your code as far as I can see. The deviation is inherent to the method.

@acer Thanks for pointing that out. I think that I was completely wrong about what the PolyhedralSets package does. I think now that it may be restricted to convex sets, but I'm not sure. Perhaps it would work if the polygon were expressed as a union of convex polygons. 

@jrive The result argument(a + b*I) = arctan(b/a) is not valid if a=0. So, if rather than assuming that your variables are real, you do

assum:= R1>0, C1>0, omega>1, vab>0;

then the end result will be Pabang = arctan(1/(omega*C1*R1)).

@C_R That's correct: simplify is executed first, and then convert. That's true whether or not the simplify is specified as a separate command. What made you think that there was a difference?

@mmcdara I don't think that a recurrence obtained using a finite number of specific coefficients can be considered a "proof". I used FormalPowerSeries instead of taylor so that I could work with arbitrary coefficients and thus get a genuine proof; specifically, (-1)^k/(2*k+1)! for sin(x) and (-4)^k/(2*k+1)! for sin(x)*cos(x). Those are the coefficients of the nonzero terms. To make them correspond to the degree d of the term, I substituted k= (d-1)/2 (note that d is always odd).

@Carl Love The 2nd part of Acer's Answer shows a better way to do the above.

@Ronan Algsubs is chaos and can't be compared to simplify with side relations which is using a well-studied bona fide algorithm called Groebner bases/Buchberger's algorithm. Like all true algorithms, it has a well-defined stopping criterion. However, its performance is very sensitive to the "monomial order" used. There are several paragraphs about this on the help page ?simplify,siderels.

@Ronan In the example that you posted, I don't understand the order (or stopping criterion) that simplify is using. Your final polynomial implicitly contains the degree-2 monomial c__2*c__3. To my mind, this should've been reduced to degree-1 by eq4. But that's not done if you don't use mindeg or call simplify again using a smaller set of equations as 2nd argument. 

@The function You're misinterpreting the plot command. Consider the command

plot([f(t), g(t), t= a..b], ...)

where the list in the 1st argument necessarily has 3 members. This command is not meant to plot f(t) and g(t) as separate functions. Rather, it plots the parametric curve {x = f(t), y = g(t), a <= t <= b}. This is a convenient way to represent curves (such as circles) that aren't properly functions.

@The function No, there is no need that the function have a real root. Consider:

eq:= (x-2)^2+(y-2)^2 = 1:
f:= solve(eq, y)[1];
plot(f, x= 1..3, scaling= constrained, view= [(0..3)$2], gridlines);

This could be rotated around the x-axis just as easily as your original.

@rcorless The integral hides a very simple conditionally convergent alternating series (shown below) akin to an alternating harmonic series. That makes a Psi representation seem quite plausible to me.

Using assuming in the above Answer wasn't necessary; option continuous handles the assumptions.

I made that derivation a bit more automatic:

restart:
#Convert F(G(g1,gk),fk) to G(F(g1,fk),gk)
SwitchOps:= (F::And(function, function &under curry(op,1)))->
local `&<`:= rcurry(op, F);
    (&<[1,0])((&<0)(&<[1,1], &<(2..)), &<[1, 2..])
:
J:= n-> Int(1/(1+x^n), x= 0..1, continuous):
S:= thaw(applyop(convert, 1, subs(x^n= freeze(x^n), J(n)), FPS));
                    /  infinity              
                   | 1  -----                
                   |     \               k   
                   |      )        k / n\    
                   |     /     (-1)  \x /  dx
                   |    -----                
                  /     k = 0                
                   0                         

S2:= SwitchOps(S);
                 infinity                     
                  -----  /  /1               \
                   \     | |             k   |
                    )    | |       k / n\    |
                   /     | |   (-1)  \x /  dx|
                  -----  \/                  /
                  k = 0    0                  

S3:= applyop(value, 1, S2);
                        infinity       
                         -----         
                          \          k 
                           )     (-1)  
                          /     -------
                         -----  k n + 1
                         k = 0         

F:= unapply(simplify(value(S3)), n);
      n-> (Psi(1/2/n + 1/2) - Psi(1/2/n))/2/n

#Verification (only valid for n > 0):
evalf(F(9) = J(9));
                  0.9320304240 = 0.9320304242

value(limit(F(n), n= 0, right) = J(0));
                             1   1
                             - = -
                             2   2

@hieudeptrai Enter ?DEplot at a command prompt to get to the help page.

@pazduha Then do

xk = eval(xk, solve({eq||(1..3)}, {xk, xdot, xdot2}));