Laplace's equation with mixed boundary condition using separation of variables

The primary goal is to solve the pde symbolically.

Then forget separation of variables, try finiteFourierSinTransform. First interpret the PDE and corresponding b.c.s to Mathematica code.

With[{T = T[x, y]}, eq = D[T, x, x] + D[T, y, y] == 0;
  bcx = {T == T1 /. x -> 0, T == T2 /. x -> a};
  bcy = {T == T2 /. y -> 0, k D[T, y] == h (T - T3) /. y -> b}];

Then introduce the transform $T(x,y)=u(x,y)+frac{T_2-T_1}{a}x+T_1$ to make the b.c.s in $x$ direction homogeneous. This isn't actually necessary, but will make the resulting series solution converge faster.

transform = T -> Function[{x, y}, u[x, y] + (T2 - T1)/a x + T1];

{neweq, newbcx, newbcy} = {eq, bcx, bcy} /. transform // Simplify

Eliminate derivative in $x$ direction using finiteFourierSinTransform:

tset = finiteFourierSinTransform[{neweq, newbcy}, {x, 0, a}, n] /. Rule @@@ newbcx

Solve the resulting ODE system:

tsol = u[x, y] /. 
   First@DSolve[tset /. HoldPattern@finiteFourierSinTransform[a_, __] :> a, u[x, y], y] //
   Simplify

Transform back:

sol = T[x, y] /. transform /. 
  u[x, y] -> inverseFiniteFourierSinTransform[tsol, n, {x, 0, a}]

Notice I've used C to denote $infty$ in sol.

To check its validity, compare it to the numeric solution:

xR = 13; yR = 17;

para = {k -> 2, h -> 3, T1 -> 5, T2 -> 7, T3 -> 11, a -> xR, b -> yR};

test[x_, y_] = sol /. para /. C -> 10 // ReleaseHold;

nsol = NDSolveValue[{eq[[1]] == NeumannValue[-h/k (T[x, y] - T3), y == b], bcx, 
    bcy[[1]]} /. para, T, {x, 0, xR}, {y, 0, yR}]

Manipulate[Plot[{nsol[x, y], test[x, y]}, {x, 0, xR}, 
  PlotStyle -> {Automatic, {Red, Dashed}}, PlotLegends -> {"FEM", "Series(10 terms)"}, 
  PlotRange -> {0, 15}], {y, 0, yR}]

Sorry, I found minor issue actually after I verified with Maple. When breaking the problem into 4 problems, to make life easier for DSolve, the top BC, should be k*Derivative[0, 1][T][x, b] == h*(T[x, b] - T3) when non homogeneous, but for the homogeneous case, it should be k*Derivative[0, 1][T][x, b] == h*(T[x, b]) and not k*Derivative[0, 1][T][x, b] == 0 as I had it below. i.e. only T3 should be set to zero.

I noticed this when I compared Maple's solution to Mathematica's and found very small difference in numerical values.

But now DSolve now can't solve it any more :(

So will post Maple's solution below, and post the corrected Mathematica solution, which breaks the problems into 4, but now it does not solve it. So may be you can try NDSolve in Mathematica for this. May be in V 12.2 DSolve can do it.

Maple solution

restart;
pde := diff(T(x, y),x$2)+diff(T(x,y),y$2)=0;
bc := T(0, y) = T1, T(a, y) = T2, T(x, 0) =T2, k*D[2](T)(x,b)= h*(T(x, b) - T3);
sol1:=simplify(pdsolve([pde, bc], T(x, y)) assuming a>0,b>0);

which gives

T(x,y) = (Sum(-2*(-(-T1+T2)*(-Pi*k*n+a*h)*exp(n*Pi*(2*b-y)/a)+a*((T2-T3)*(-1)^n
-T1+T3)*h*exp(n*Pi*(b-y)/a)-a*((T2-T3)*(-1)^n-T1+T3)*h*exp(n*Pi*(y+b)/a)+exp(n/
a*Pi*y)*(-T1+T2)*(Pi*k*n+a*h))*sin(n/a*Pi*x)/Pi/n/((-Pi*k*n+a*h)*exp(2*n/a*Pi*b
)-Pi*n*k-a*h),n = 1 .. infinity)*a+T1*a+x*(-T1+T2))/a

Corrected Mathematica solution

ClearAll[T, x, y, T1, T2, T3, k, h];
pde = Laplacian[T[x, y], {x, y}] == 0;
bc = {T[0, y] == 0, T[a, y] == 0, T[x, 0] == 0, k*Derivative[0, 1][T][x, b] - h*T[x, b] == -h*T3}
sol1 = DSolve[{pde, bc}, T[x, y], {x, y}, Assumptions -> {a > 0, b > 0}]

bc = {T[0, y] == 0, T[a, y] == 0, T[x, 0] == T2, k*Derivative[0, 1][T][x, b] - h*T[x, b] == 0}
sol2 = DSolve[{pde, bc}, T[x, y], {x, y},   Assumptions -> {a > 0, b > 0}]

bc = {T[0, y] == 0, T[a, y] == T2, T[x, 0] == 0,   k*Derivative[0, 1][T][x, b] - h*T[x, b] == 0}
sol3 = DSolve[{pde, bc}, T[x, y], {x, y},   Assumptions -> {a > 0, b > 0}]

(* no solution *)

bc = {T[0, y] == T1, T[a, y] == 0, T[x, 0] == 0, k*Derivative[0, 1][T][x, b] - h*T[x, b] == 0}
sol4 = DSolve[{pde, bc}, T[x, y], {x, y},   Assumptions -> {a > 0, b > 0}]

(* no solution *)

So can't add them, since 2 solutions could not be found. need all 4.

Original answer below (but contain error in BC)

You can solve it as follows.

DSolve can not handle more than one edge in Laplacian in 2D being non homogeneous at same time.

Due to linearity, the problem is broken into 4 problems, were one edge is non homogeneous at time. Then the 4 solutions are added.

So instead of doing

ClearAll[T, x, y, T1, T2, T3, k, h];
pde = Laplacian[T[x, y], {x, y}] == 0;
bc = {T[0, y] == T1, T[a, y] == T2, T[x, 0] == T2, k*Derivative[0, 1][T][x, b] == h*(T[x, b] - T3)}
sol1 = DSolve[{pde, bc}, T[x, y], {x, y}, Assumptions -> {a > 0, b > 0}]

Where DSolve does not solve it, do the following

bc = {T[0, y] == 0, T[a, y] == 0, T[x, 0] == 0, k*Derivative[0, 1][T][x, b] == h*(T[x, b] - T3)}
sol1 = DSolve[{pde, bc}, T[x, y], {x, y},   Assumptions -> {a > 0, b > 0}];

bc = {T[0, y] == 0, T[a, y] == 0, T[x, 0] == T2,k*Derivative[0, 1][T][x, b]  == 0}
sol2 = DSolve[{pde, bc}, T[x, y], {x, y},   Assumptions -> {a > 0, b > 0}]

bc = {T[0, y] == 0, T[a, y] == T2, T[x, 0] == 0, k*Derivative[0, 1][T][x, b]  == 0}
sol3 = DSolve[{pde, bc}, T[x, y], {x, y},   Assumptions -> {a > 0, b > 0}]

bc = {T[0, y] == T1, T[a, y] == 0, T[x, 0] == 0, k*Derivative[0, 1][T][x, b] == 0}
sol4 = DSolve[{pde, bc}, T[x, y], {x, y},   Assumptions -> {a > 0, b > 0}]

And the solution is

sol = First[(T[x, y] /. sol1) 
     + (T[x, y] /. sol2) 
     + (T[x, y] /. sol3) 
     + (T[x, y] /. sol4)]