I'm not sure how to get a symbolic solution.  We can, however, symbolically solve the steady-state periodic output:

dsol_ss := dsolve({diff(x,t)+x=Et, eval(x,t=0)=eval(x,t=T)});
dsol_ss :=

        {                 2 exp(-t) (-1 + exp(-1))
        {         t - 1 + ------------------------                 t < 1
        {                       exp(-2) - 1
    x = {
        {          2 exp(-t) (-1 + exp(-1))
        { -t + 3 + ------------------------ - 2 exp(-t + 1)        1 <= t
        {                exp(-2) - 1

plot([rhs(dsol_ss),Eper],t=0..T);

Actually, solving the complete system symbolically isn't much more difficult, since Maple's dsolve returns a result for this diff-eq with a symbolic initial condition.

Sure, at least numerically:

Et := piecewise(t<=1,t,2-t):
T := 2:
Eper := eval(Et, t=t-T*trunc(t/T));
alias(x=x(t)):
deq := diff(x,t) + x = Eper:
ics := eval(x,t=0)=0:
dsol := dsolve({deq,ics},'numeric'):
plots[odeplot](dsol, [[t,x],[t,Eper]],0..5);

Are you trying to simulate (say numerically) the response of the circuit to the periodic waveform, or are you trying to solve the corresponding periodic output?  A usual method for generating a periodic numerical waveform is to use replace the t in the expression with a periodic expression in t:

eval(expr, t=t-T*trunc(t/T)): # T is the period

For example, here is how to make a periodic triangular waveform:

Et := piecewise(t<=1,t,2-t):
T := 2:
Eper := eval(Et, t=t-T*trunc(t/T)):
plot(Eper, t=0..5);

That's just the way Maple works.  It automatically distributes multiplication of numerics over sums.  You can fake the output using the empty function (``):

x := a/2+b/3:
makefrac := (``@numer)/(``@denom):
makefrac(x);
                                  (3 a + 2 b)
                                 ------------
                                      (6)
expand(%);
                                   a/2 + b/3

How 'bout

eq1 := 3*x + y - z = 3:
eq2 := x + 2*y + 4*z = -4:
sol := solve({eq1,eq2},{y,z});
                                                      13 x            5 x
                                        sol := {y = - ---- + 4/3, z = --- - 5/3}
                                                       6               6

line := subs(sol, x=t, <x,y,z>);
                                                         [     t      ]
                                                         [            ]
                                                         [  13 t      ]
                                                         [- ---- + 4/3]
                                                 line := [   6        ]
                                                         [            ]
                                                         [ 5 t        ]
                                                         [ --- - 5/3  ]
                                                         [  6         ]

subs(t=0,line);
                                                         [ 0  ]
                                                         [    ]
                                                         [4/3 ]
                                                         [    ]
                                                         [-5/3]

One could also use LinearAlgebra:-IntersectionBasis to find the coefficients:

S1 := convert(map2(coeff, lhs(eq1), [x,y,z]), 'Matrix'):
S2 := convert(map2(coeff, lhs(eq2), [x,y,z]), 'Matrix'):
N1 := NullSpace(S1):
N2 := NullSpace(S2):
IntersectionBasis([N1,N2]);
                                                         [-6 ]
                                                         [-- ]
                                                         [13 ]
                                                         [   ]
                                                        {[ 1 ]}
                                                         [   ]
                                                         [-5 ]
                                                         [-- ]
                                                         [13 ]

One possibility is to first solve the equation symbolically, before assigning values to the parameters.  In fact, it is usually best to avoid directly assigning values to parameters:

eq := K = S*F^n:
vals := [K=2, S=3, n=-2.4]:
sol := solve(eq, F);
evalf(eval(sol, vals));

Another approach is to use the RealDomain package to restrict solutions to reals:

RealDomain:-solve(eval(eq,vals),F);

It isn't clear to me that acer's suggestion does what you ask.  That is, because simplify(u*Dirac(u)) returns 0, applying simplify to den returns 0, while you asked for v*Dirac(v).  To achieve that, you might use frontend and specify that the function Dirac(u) not be frozen:

den := a*u*Dirac(u)+b*u*u*Dirac(u)+c*u^(3/2)*Dirac(u)+v*Dirac(v);
                              2               (3/2)
     den := a u Dirac(u) + b u  Dirac(u) + c u      Dirac(u) + v Dirac(v)

frontend(simplify, [den], [{`*`,`+`,`^`},{Dirac(u)}]);
                                  v Dirac(v)

This might be appliable to your followup question.

Good question.  It is a habit from using the older matrix structure.  There eval was needed because matrices (small-m), as well as tables, arrays, procedures, and modules (records) use last-name-evaluation, which means that if M were a (small-m) matrix, then returning M would actually return the local name M rather than the matrix M.  Because rtables (of which Matrix is a specific type) do not use last-name-evaluation this isn't necessary. 

See the help page ?last_name_eval for details.

FYI, your procedure had a couple of flaws.  First, it never assigned a Matrix to M.  Second, it never returned the computed Matrix.  Here is a corrected version.  Of course, the methods that have been suggested are more convenient.

mat:= proc(a,b,c,N)
local i,j,M:
    M := Matrix(N,N);
    for i from 1 to N do
        for j from 1 to N do
            if j>(i+1) then M[i,j]:= a:
            elif i>(j+1) then M[i,j]:= c:
            else M[i,j]:= b:
            end if:
        end do:
    end do:
    eval(M);
end proc:
mat(a,b,c,4);

> Matrix(4,4,(i,j)->`if`(j>i+1,a,`if`(i>j+1,c,b)));   
                                  [b    b    a    a]
                                  [                ]
                                  [b    b    b    a]
                                  [                ]
                                  [c    b    b    b]
                                  [                ]
                                  [c    c    b    b]

eq1 := x^3 + y^3 = 3*x^2*y:
eq2 := subs(y=y(x),eq1);
                                 3       3      2
                         eq2 := x  + y(x)  = 3 x  y(x)

isolate(diff(eq2,x),diff(y(x),x));
                                        2
                          d         -3 x  + 6 x y(x)
                          -- y(x) = ----------------
                          dx               2      2
                                     3 y(x)  - 3 x

subs(y(x)=y,rhs(%));
                                     2
                                 -3 x  + 6 x y
                                 -------------
                                     2      2
                                  3 y  - 3 x

alias(y=y(x)):
eq := x^3 + y^3 = 3*x^2*y:
isolate(diff(eq,x),diff(y,x));

                                        2
                             d      -3 x  + 6 x y
                             -- y = -------------
                             dx         2      2
                                     3 y  - 3 x

First off, your definition of coprime is nonstandard (putting it politely).

I'd try

iscoprime := (k,j) -> evalb(igcd(k,j)=1):
L := [seq(1..65)]:
select(iscoprime, L, 66):

If you suspect that the answer is a rational polynomial, you can try to solve for the coefficients

f := 1/(k^2+3*k+1):

nn := 1:
nd := 3:

g := add(a||i*k^i,i=0..nn)/add(b||i*k^i,i=0..nd);
                                     a0 + a1 k
                        g := -------------------------
                                             2       3
                             b0 + b1 k + b2 k  + b3 k


eqs := {seq(numer(eval(g-f,k=kk)),kk=1..nn+nd+2)}:
sol := solve(eqs);
sol := {a0 = -3 b3 + b2, a1 = b3, b0 = -3 b3 + b2, b1 = -8 b3 + 3 b2, b2 = b2,

    b3 = b3}

normal(subs(sol,g));
                                      1
                                 ------------
                                  2
                                 k  + 3 k + 1

Here's a useful heuristic.  Given a problem that you cannot solve, create an easier version of it that you can, than attempt to extend the solution.  That is the reason you were asked to solve for p(2) and p(3) before p(n).  However, you can make it even simpler.  The general problem is, given N days in a year ... find p(n), where n is the number of  people ....   Suppose N=1.  What is p(2)?   p(3)?  p(n)?  Now suppose N=2.

That the ring and the sphere have different Gaussian curvatures directly proves that there will be a "bucklement", however, if you replace the sphere with a cone you will also get a "bucklement" even though the cone has zero Gaussian curvature.  You can, however, wrap the ring around the cone (rather than encircling it and then trying to squash it flat).