How to cleverly solve this kind of equations? 
Suppose we have the following equations where $0<p<1$:
\begin{align} a_{i,j}&=pb_{i+1,j} + (1-p)b_{0,j} \\
 b_{i,j}&=pa_{i,j+1}+(1-p)a_{i,0} \end{align} with the boundary
  conditions: \begin{align} b_{3,j}=1 \qquad\text{  where $0\le i< 3$}\\
 a_{i,3}=0 \qquad\text{  where $0\le j< 3$} \end{align}
  We are asked to get $a_{0,0}$

I can solve it by getting a set of equations such as:
\begin{align}
a_{0,0} & = pb_{1,0} + (1-p) b_{0,0}\\
b_{1,0} & = pa_{1,1} + (1-p) a_{1,0}\\
b_{0,0} & = pa_{0,1} + (1-p) a_{0,0}\\
a_{1,1} & = pb_{2,1} + (1-p) b_{0,1}\\
\vdots
\end{align}
however this seems too complicated and very easy to calculate wrong. If we let the number $3$ in the boundary condition to be a much larger number $n$, then this method is impossible to solve by hand. 
I want to know if there is a better way or more analytical way to solve this kind of equation?
 A: A general method that works for any size $n$ and is not prone to errors is to use a computer algebra system (such as Mathematica) to directly solve the system for any size $n$. Given that this is a well-posed finite linear system of equations, I don't see any clever tricks; if you want to solve such a system, you (or the computer) has to do the tedious work of solving it using elimination of variables until you get the result. Whether you call this Gaussian elimination or LU decomposition, you have to substitute the equations into each other and keep track of your operations without error in order to arrive at the correct result. In much more complicated problems (such as solving PDEs such as the wave or heat equation or Maxwell's equations), once you reduce things to a finite system of linear equations, you let the computer solvers take over. The fact that your problem is of this simple class to begin with leads me to believe there really aren't any clever tricks to be played here.
For the approach, I transcribed the two equations you presented into Mathematica, then set it up to explicitly calculate all the equations for $i,j\in\{0,1\ldots n\}$, enforced the boundary conditions for $i=n+1$ and $j=n+1$, then told it to solve the equations. For $n=2$, the solution is:
$$
a_{0,0} =  -\frac{p^8-2 p^7+p^6-p^4+p^3-1}{p^9-4 p^8+6 p^7-4 p^6+3 p^4-3 p^3+2}
$$
A happy consequence is that all the unknowns are solved for simultaneously by Mathematica's algorithms.
Here is a plot of $a_{0,0}$ as a function of $p$:

And here is the Mathematica code I used:
(*setup the right hand side functions*)
a[i_,j_] := p Subscript[b, i+1,j] + (1-p) Subscript[b, 0,j]
b[i_,j_]:=p Subscript[a, i,j+1] + (1-p) Subscript[a, i,0]
n=2; (*works for any order n*)
(*set the variables equal to the right hand sides, enforce the boundary conditions, and collect into one list*)
aEqs = Table[Subscript[a, i,j]==a[i,j],{i,0,n},{j,0,n}]/.{Subscript[a, i_,n+1]->0,Subscript[b, n+1,j_]->1}//Flatten[#,1]&;
bEqs = Table[Subscript[b, i,j]==b[i,j],{i,0,n},{j,0,n}]/.{Subscript[a, i_,n+1]->0,Subscript[b, n+1,j_]->1}//Flatten[#,1]&;
Eqs = Join[aEqs,bEqs];
(*collect all the variables into one list*)
aVars = Table[Subscript[a, i,j],{i,0,n},{j,0,n}]//Flatten[#,1]&;
bVars = Table[Subscript[b, i,j],{i,0,n},{j,0,n}]//Flatten[#,1]&;
Vars = Join[aVars,bVars];
(*Run the solver*)
Soln = Subscript[a, 0,0]/.Solve[Eqs,Vars][[1]]//Simplify;
TraditionalForm[Soln]
Plot[Soln,{p,0,1}]

