$(\partial_{tt}-\nabla^2+\partial_t)f=g,\quad (\partial_t-\nabla^2+b)g=\partial_t f$ Hi I am looking for complete solutions for $f(r,t),g(r,t)$ given in the coupled linear partial differential equations below:
$$
(\partial_{tt}-a\nabla^2+b\partial_t)f(r,t)=bg(r,t)
$$
$$
(\partial_t-c \nabla^2+b)g(r,t)=b\partial_t f(r,t)
$$
Initial and boundary conditions are given by
$$
\partial_t f(0,t)=0,\ \partial_t f(R,t)=d\cos(\omega t)
$$
$$
g(0,t)=0, \  g(R,t)=d\cos (\omega t)
$$
where $a,b,c,d,R,\omega>0.$  Note
$$
\nabla^2\equiv \frac{1}{r}\partial_r(r\partial_r)-\frac{1}{r^2}=\partial_{rr}+\frac{1}{r}\partial_r -\frac{1}{r^2}.
$$
Thank you!  Some comments:
If the right hand side of both equations are zero, (equations become homogenous), then the kernel of both linear operators are known and are in terms of Bessel functions $J_1$ (we're in a cylindrical geometry hence the Laplacian like term) times oscillating functions of time: $\sin \omega t,\cos \omega t$. 
 A: This is sketch of how to solve your problem.  I still don't understand the boundary conditions/ initial conditions as stated.  I have set $a,b,d$ to one, and have made generic boundary conditions.
Let $\Omega \subset \mathbb{R}^d$ be a bounded open domain.  We are interested in solving the following coupled systems of PDE.  This is pretty hairy.
$$
\begin{align*}
0 &= (\partial_{tt} - \nabla^2 + \partial_t) f - g&\quad \mathrm{for}\  (x,t)\in \Omega \times \mathbb{R}^+ \\
0 &= (\partial_t-\nabla^2 + 1) g - \partial_t f&\quad \mathrm{for} (x,t) \in \Omega \times \mathbb{R}^+ \\
f &= 0 &\quad  \mathrm{for} (x,t)\in \Omega\times \{0\}\\ 
f_t &= 0&\quad \mathrm{for} (x,t)\in \Omega\times \{0\}\\
g &= 0& \quad \mathrm{for} (x,t)\in \Omega\times \{0\}\\
f_t &= s_1 &\quad \mathrm{for} (x,t)\in \partial\Omega \times \mathbb{R}^+\\
g &= s_2 &\quad \mathrm{for} (x,t)\in \partial\Omega \times \mathbb{R}^+
\end{align*}
$$
The fist two equations are the PDE.  The second three are initial conditions, and the last two are boundary conditions.  
Step #1 Homogenize the PDE's boundary conditions.  We need to turn the boundary conditions into "forcing terms"
Let $f_{bc}$ be the solution to the following PDE
\begin{align*}
-\nabla^2 f_{bc} &= 0 &\quad \mathrm{for} (x,t)\in \Omega \times \mathbb{R}^+\\
\frac{\mathrm{d}}{\mathrm{d}t} f_{bc} &= s_1 & \quad \mathrm{for} (x,t)\in \partial\Omega\times \mathbb{R}^+
\end{align*}
With simple boundary conditions and simple geometry this $f_{bc}$ might just be a constant in space and vary in time.
Let $g_{bc}$ be the solution of a similar PDE
\begin{align*}
-\nabla^2 g_{bc} &= 0&\quad \mathrm{for} (x,t)\in \Omega \times \mathbb{R}^+\\
g_{bc} &= s_2& \quad \mathrm{for} (x,t)\in \partial\Omega \times \mathbb{R}^+.
\end{align*}
Now we have a set of new PDE's.
$$
\begin{align*}
-(\partial_{tt} + \partial_t )f_{bc} + g_{bc}&= (\partial_{tt} - \nabla^2 + \partial_t) (f-f_{bc}) -(g-g_{bc})&\quad \mathrm{for}\  (x,t)\in \Omega \times \mathbb{R}^+ \\
-(\partial_t +1 )g_{bc} + (f_{bc})_t&= (\partial_t-\nabla^2 + 1) (g-g_{bc}) - \partial_t (f-f_{bc})&\quad \mathrm{for} (x,t) \in \Omega \times \mathbb{R}^+ \\
f-f_{bc} &= -f_{bc} &\quad  \mathrm{for} (x,t)\in \Omega\times \{0\}\\ 
f_t-(f_{bc})_t &= -(f_{bc})_t&\quad \mathrm{for} (x,t)\in \Omega\times \{0\}\\ 
g-g_{bc} &= -g_{bc}& \quad \mathrm{for} (x,t)\in \Omega\times \{0\}\\
f_t-(f_{bc})_t &= 0 &\quad \mathrm{for} (x,t)\in \partial\Omega \times \mathbb{R}^+\\
g-g_{bc} &= 0 &\quad \mathrm{for} (x,t)\in \partial\Omega \times \mathbb{R}^+
\end{align*}
$$
Notice what we have done, we have changed a PDE with with inhomogeneous boundary conditions and no "forcing functions" into one with homogeneous boundary conditions with forcing functions.
Step #2.
Denote $\{\Phi_i\}_{i=1}^\infty$ as a sequence of eigenfunctions for the associated homogeneous Poisson equation, ordered by increasing eigenvalue.  We mean to say that $\int_{\Omega} \Phi_i \Phi_j \mathrm{d} x = \delta_{i,j}$ and $\int_{\Omega} (\nabla \Phi_i) \cdot (\nabla \Phi_j) \mathrm{d}x = \delta_{i,j} \lambda_i^2$.  (For the cylindrical case we have $\Phi_i$ are the Bessel functions.)
Now, these eigenfunctions are complete so we can expand $(f-f_{bc})$ as
$$
(f-f_{bc})(x,t) = \sum_{i=1}^\infty \tilde{f}_i(t) \Phi_i(x)
$$
and $$(g-g_bc)(x,t) = \sum_{i=1}^\infty \tilde{g}_i(t) \Phi_i(x)$$
Now multiplying the PDE by $\Phi_j$ integrating by parts the Laplacian term we get
$$
\underbrace{\int_{\Omega} \Phi_i ((-\partial_{tt} - \partial_t )f_{bc} + g_{bc} )\, \mathrm{d} x}_{\alpha_i(t)} = (\tilde{f}_i)_{tt} + (\tilde{f}_i)_t + \lambda_i^2 \tilde{f}_i - \tilde{g}_i
$$
and
$$
\underbrace{\int_{\Omega} \Phi_i ((-\partial_t+1) g_{bc} +  \partial_t f_{bc})\mathrm{d}x }_{\beta_i(t)} = (\tilde{g}_i)_t+\lambda_i^2 \tilde{g}_i + \tilde{g}_i - (\tilde{f}_i)_t
$$
The beauty is that, at this point $\alpha_i$ and $\beta_i$ are all known, and  $\tilde{f}_i$ only couples with $\tilde{g}_i$.
Step #3.  Solve ODE.
We write a state-space representation of the ODE's.
$$
\begin{align*}
\frac{\mathrm{d}}{\mathrm{d} t}\left\lbrack
\begin{array}{c}
\tilde{f}_i\\
(\tilde{f}_i)_t\\
\tilde{g}_i 
\end{array}\right\rbrack + \left\lbrack
\begin{array}{ccc}
0&-1&0\\
\lambda_i^2& 1& -1\\
0& -1 & \lambda_i^2+1
\end{array}
\right\rbrack
\left\lbrack
\begin{array}{c}
\tilde{f}_i\\
(\tilde{f}_i)_t\\
\tilde{g}_i
\end{array}
\right\rbrack
= \left\lbrack
\begin{array}{c}
0\\
\alpha_i\\
\beta_i\\
\end{array}
\right\rbrack
\end{align*}
$$
Solve this system how ever you want.  This solves $f-f_{bc}$ and $g-g_{bc}$ so to get $f$ and $g$ you just add the boundary functions.  Do you see how to handle the initial conditions?  I can explain any steps that I skipped or ones that are unclear.
A: This is an extension of the previous answers :
$(\partial_{tt}+\partial_t-\nabla^2)f(r,t)=0$
PDE $(\partial_{tt}+a\partial_t-b\nabla^2)f(r,t)=0$
In the calculus below, some results which have already been obtained are re-used without more explanation.
 

